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ABSTRACT 

We present a method for modelling star-forming clouds that combines two different models 
of the thermal evolution of the interstellar medium (ISM). In the combined model, where the 
densities are low enough that at least some part of the spectrum is optically thin, a model 
of the thermodynamics of the diffuse ISM is more significant in setting the temperatures. 
Where the densities are high enough to be optically thick across the spectrum, a model of 
fiux limited diffusion is more appropriate. Previous methods either model the low-density 
interstellar medium and ignore the thermal behaviour at high densities (e.g. inside collapsing 
molecular cloud cores), or model the thermal behaviour near protostars but assume a fixed 
background temperature (e.g. ^ 10 K) on large-scales. Our new method treats both regimes. 
It also captures the different thermal evolution of the gas, dust, and radiation separately. We 
compare our results with those from the literature, and investigate the dependence of the 
thermal behaviour of the gas on the various model parameters. This new method should allow 
us to model the ISM across a wide range of densities and, thus, develop a more complete and 
consistent understanding of the role of thermodynamics in the star formation process. 

Key words: astrochemistry - hydrodynamics - ISM: general - methods: numerical - radia¬ 
tive transfer - stars: formation. 


1 INTRODUCTION 

The thermal behaviour of interstellar gas is of fundamental impor¬ 
tance for star formation. To initiate star formation, the thermal pres¬ 
sure must be insufficient to support the gas against gravitational 
collapse (Jeans 1929). Further evolution also depends on the ther¬ 
modynamics and density structure, with a variety of different out¬ 
comes being possible such as the formation of a single polytrope or 
fragmentation into many clouds (Hoyle 1953; Hunter 1962; Layzer 
1963; Tohline 1980; Rozyczka 1983). Larson (1985, 2005) empha¬ 
sised the importance of the relationship between temperature and 
density in the interstellar medium, and Whitworth, Boffin & Fran¬ 
cis (1998) emphasised the importance of the transition from molec¬ 
ular cooling to dust cooling. The typical decrease in gas tempera¬ 
ture from Tg ~ 1000 K at number densities of hydrogen nuclei of 
nn ~ 1 cm“^ to Tg ~ 10 K at densities of nn ~ 10^ cm“^ pro¬ 
motes fragmentation, while the transition to an isothermal regime 
or temperature that increases with density above nn ~ 10^ cm“^ 
may help to produce a characteristic stellar mass (Larson 1985, 
2005; Jappsen, Klessen, Larson, Li & Mac Low 2005; Bonnell, 
Clarke & Bate 2006; Elmegreen, Klessen & Wilson 2008). Chemi¬ 
cal reactions that control the abundances of gas phase coolants, and 
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therefore radiative equilibrium, may affect this transition and the 
formation of molecular cloud cores. At even higher densities, deep 
within collapsing molecular cloud cores, different phases of col¬ 
lapse are believed to occur as the gas becomes optically thick and 
later as molecular hydrogen dissociates (Larson 1969). The treat¬ 
ment of radiative transfer can be crucial for determining whether 
the core fragments into multiple protostars or not (Boss et al. 2000; 
Krumholz 2006; Whitehouse & Bate 2006). 

Previous radiation hydrodynamical simulations of star cluster 
formation demonstrate the importance of radiative feedback from 
protostars to correctly capture fragmentation and produce realistic 
numbers of brown dwarfs (Bate 2009; Offner et al. 2009). Whereas 
barotropic calculations of star cluster formation produce charac¬ 
teristic stellar masses that depend on the initial Jeans mass in the 
cloud (Bate & Bonnell 2005), Bate (2009) showed that radiative 
feedback reduces this dependence on the initial conditions. This 
may help to explain the apparent universality of the stellar initial 
mass function (IMF), at least in recent epochs (see also Krumholz 
2011). Indeed, some recent radiation hydrodynamical simulations 
of star cluster formation have been quite successful in reproduc¬ 
ing the observed IMF and other properties of stellar systems (Bate 
2012, 2014; Krumholz, Klein & McKee 2012). 

However, to date, radiation hydrodynamical simulations of 
star formation have been restricted to modelling dense molecu- 
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lar clouds and assuming that the only significant sources of ra¬ 
diation are the protostars themselves. In this case, the gas tem¬ 
perature at large distances from the protostars is assumed to be 
10 K. However, this approximation is generally invalid at densi¬ 
ties nn ^ 10^ cm“^ with solar metallicities, and is invalid at even 
higher densities at lower metallicities, if at all. 

Informed by the vast literature on the physics, thermodynam¬ 
ics, and chemistry of the diffuse interstellar medium (ISM) (see 
the thorough review provided by Tielens 2005) and collapsing 
clouds with different metallicities (e.g. Omukai 2000; Omukai et al. 
2005), interest has grown in studying the formation and evolution 
of molecular clouds using three-dimensional hydrodynamical cal¬ 
culations coupled with chemistry. Initial caleulations treated only 
the formation of molecular hydrogen (Dobbs, Bonnell & Pringle 
2006; Dobbs & Bonnell 2007; Glover & Mae Low 2007a,b; Micic 
et al. 2012). More recent work has included more complex chem¬ 
istry, in particular that involving carbon and oxygen (Glover et al. 
2010; Glover & Mac Low 2011; Glover & Clark 2012b,a,e; Clark 
et al. 2012). Along with the formation mechanisms of molecular 
clouds and their chemical and temperature distributions, some of 
these calculations have begun to investigate star formation, using 
sink particles to replace collapsing regions of gas (in particular the 
models of Glover & Clark). At the same time, there have been sev¬ 
eral studies of the chemical and thermal evolution of low metallic- 
ity clouds and their star formation (e.g. Dopcke et al. 2011, 2013; 
Glover & Clark 2012c). 

Therefore, at the present time, the star formation community 
has two different classes of hydrodynamieal models for star forma¬ 
tion. One class treats the low-density thermochemical evolution in 
some detail, but does not treat radiative feedback from protostars. 
The other the ignores the complieated physics of the diffuse ISM, 
which is particularly important at low densities and low metallic¬ 
ities, but includes the radiative effects of protostars. The goal of 
this paper is to make an attempt to bridge this gap, by combining a 
model for the physics of the low-density interstellar medium with 
a method for modelling radiative transfer. Our main purpose is to 
develop a model that does a reasonable job of modelling the ther¬ 
modynamics at both low and high densities in star-forming clouds, 
with metallicities as low as 1/100 solar. The aim is not to produce 
a detailed chemical model. Rather, we wish to implement the sim¬ 
plest possible chemical model that will provide the abundances of 
the major coolants of the gas that are necessary to calculate realistic 
temperatures. 

It turns out that a method very similar to that which we present 
here has recently been developed by Pavlyuchenkov & Zhilkin 
(2013) and Pavlyuchenkov et al. (2015). Although we developed 
our method independently, each of the methods is based on extend¬ 
ing a method of radiative transfer that is valid at high densities to 
include effeets that are important at low densities. Eaeh of the meth¬ 
ods treat gas and dust temperatures separately. Pavlyuchenkov & 
Zhilkin (2013) base their method on the diffusion approximation 
for radiative transfer and add heating and eooling terms relevant 
at lower densities to model the collapse of dense molecular cloud 
cores. Pavlyuchenkov et al. (2015) replace the diffusion approxi¬ 
mation with fiux-limited diffusion like we use in this paper. They 
do not include cooling processes relevant at very low-densities that 
we include (e.g. electron recombination and fine structure emission 
from atomic oxygen), and they only perform one-dimensional cal¬ 
culations, but the underlying methods are very similar. 

This paper is primarily a method paper where we describe our 
implementation and demonstrate its performance in a wide variety 
of tests, comparing our results to those of others who have per¬ 


formed similar caleulations (sometimes using mueh more eomplete 
and/or complex models). However, we also explore the effects of 
varying many of the parameters that go into the model in order to 
better understand which physical processes may be most important 
for affecting star formation. Large-scale star formation calculations 
are beyond the scope of this paper, but we hope to use this new 
method to perform such calculations in the future. 

In Section 2, we provide the fundamental equations and as¬ 
sumptions that go into our model. The implementation of these 
equations into a smoothed partiele hydrodynamics (SPH) code is 
described in Section 3. We present the results from our test calcula¬ 
tions, and compare our results with those in the literature in Section 
4. Finally, we draw our conclusions in Section 5. 


2 METHOD 

2.1 The flux-limited diffusion approximation 

In a frame co-moving with the fluid, and assuming local thermal 
equilibrium (LTE), the equations governing the time-evolution of 
radiation hydrodynamics (RHD) can be written 

g^+pv-v = 0, (1) 



-Vp+ , 

C 


( 2 ) 


~ ^ ~ VvrP + AiTKppB — ck,epE , (3) 

P^ (p) ^ ^ ~ ^TTHippB + CHiEpE , (4) 

(Mihalas & Mihalas 1984; Turner & Stone 2001; Whitehouse, Bate 
& Monaghan 2005). In these equations, D/Dt = d/dt v ■ V 
is the convective derivative. The symbols p, e, v and p repre¬ 
sent the material mass density, energy density, veloeity, and scalar 
isotropic pressure, respectively, and c is the speed of light. The total 
frequency-integrated radiation density, momentum density (flux) 
and pressure tensor are represented by E, F, and P, respectively. 
The assumption of LTE allows the rate of emission of radiation 
from the matter in equations 3 and 4 to be written as the frequency- 
integrated Planck function, B. Equations 2-4 have been integrated 
over frequency, leading to the flux mean total opacity xf, and the 
Planck mean and energy mean absorption opacities, kp and Kp. 
The total opacity, x. is the sum of components due to absorption k, 
and scattering a. 

Taking an ideal equation of state for the gas pressure p = 
(7 — l)up, where u = e/p is the specific energy of the gas, the 
temperature of the gas is Tg = (7 — where p is the 

mean molecular weight of the gas and IZ is the gas constant. The 
frequency-integrated Planck function is given hy B — (crs/TrjTg, 
where ctb is the Stefan-Boltzmann constant. The radiation energy 
density also has an associated temperature Tr from the equation 
E = 4(7BTpc. 

A common approximation to make in radiation hydrodynam¬ 
ics is the so-called flux-limited diffusion approximation. For an 
isotropie radiation field P = FI/3. The Eddington approximation 
assumes this relation holds everywhere and implies that, in a steady 
state 

F=-^VE. (5) 
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This expression gives the correct flux in optically thick regions, 
where XP is large. However in optically thin regions where xp ^ 0 
the flux tends to inflnity whereas in reality 1-^1 ^ cE. Flux-limited 
diffusion solves this problem by limiting the flux in optically thin 
environments to always obey this inequality. Levermore & Pom- 
raning (1981) wrote the radiation flux in the form of Pick’s law of 
diffusion as 


F = -DVE, (6) 

with a diffusion constant given by 

D=—. (7) 

XP 

The dimensionless function \{E) is called the flux limiter. The ra¬ 
diation pressure tensor may then be written in terms of the radiation 
energy density as 

P= iE, (8) 


where the components of the Eddington tensor, f, are given by 

f= i(l-/)I+^(3/-l)nn, (9) 

where fi — VF^/1| is the unit vector in the direction of the radi¬ 
ation energy density gradient and the dimensionless scalar function 
f{E) is called the Eddington factor. The flux limiter and the Ed¬ 
dington factor are related by 

/ = A + A^i^^ (10) 


where R is the dimensionless quantity R = \VE\/(xpE)- 

Equations 6 to 10 close the equations of RHD. However, we 
must still choose an expression for the flux limiter, A. In this paper, 
we choose the flux limiter of Levermore & Pomraning (1981) 


\{R) 


2-\-R 

6 + 3J7 + • 


( 11 ) 


There are a number of problems with using the flux-limited 
diffusion approximation to model radiative processes in star for¬ 
mation. Many of these are a direct result of the approximations 
made in deriving the method. The diffusion approximation is good 
at high optical depths, but the directional behaviour of anisotropic 
radiation at low optical depths is modelled poorly. This means, for 
example, that shadowing is not reproduced. The flux limiter gives 
the correct limiting behaviour for the propagation rate in the opti¬ 
cally thin and thick regimes, but at intermediate optical depths is 
dependent on the arbitrary choice of the flux limiter. Einally, most 
implementations take the grey approach, ignoring the frequency de¬ 
pendence of the radiative transfer in favour of mean opacities (as 
discussed above). This is likely to be a good approximation when 
most of the energy is in long-wavelength radiation, but is a poor 
approximation near hot massive stars (Wolflre & Cassinelli 1987; 
Preibisch, Sonnhalter & Yorke 1995; Yorke & Sonnhalter 2002; 
Edgar & Clarke 2003). 

Despite these drawbacks, flux-limited diffusion is expected 
to do a reasonable job of modelling radiative transfer in dense 
star-forming regions that produce low-mass stars (e.g. Bate 2009; 
Offner et al. 2009; Bate 2012) and it is computationally efficient. 


2.2 The diffuse interstellar medium 

However, other problems arise when modelling low-density envi¬ 
ronments. Eor a start, the above equations treat radiation and mat¬ 
ter, but do not distinguish between different types of matter which 


may have different temperatures. In particular, at solar metallicities, 
dust and gas temperatures are only well coupled (by collisions) at 
gas densities above 10^ cm“^ (Burke & Hollenbach 1983; Gold¬ 
smith 2001; Glover & Clark 2012b). At lower densities, the gas and 
dust temperatures can be very different from one another with the 
gas temperature typically exceeding that of the dust. There are also 
sources of heating that affect the matter other than work done on the 
gas and the radiative interaction between the matter and the radia¬ 
tion held that are included in equation 4. Cosmic rays heat the gas 
by direct collision (Goldsmith, Habing & Eield 1969; O’Donnell & 
Watson 1974; Black & Dalgarno 1977; Goldsmith & Danger 1978). 
Ultraviolet photons heat the gas indirectly through the photoelec¬ 
tric release of hot electrons from dust grains (Draine 1978; Bakes & 
Tielens 1994). Absorption and emission of radiation that is highly 
frequency dependent is also problematic in a grey treatment. In 
particular, gas cooling occurs by atomic and molecular line cool¬ 
ing which may be optically thin due to Doppler shifts even when 
static clouds would be optically thick. The external interstellar radi¬ 
ation (ISR) held (Habing 1968; Witt & Johnson 1973; Draine 1978; 
Black 1994), attenuated by dust extinction, also heats the dust. 


2.3 Combining a diffuse interstellar medium model with 
flux-limited diffusion 

In this section, we present a method to combine a model of the 
radiative equilibrium of the diffuse ISM with flux-limited diffusion 
to determine the gas, dust, and radiation temperatures in both low- 
density and high-density regions of molecular clouds. We include 
treatments for all of the effects mentioned in the previous section. 
We also allow for heating due to H 2 formation. We do not treat 
photoionisation or heating from X-rays. 

We use equations 2-4 to model the continuum radiation held 
and the gas, and add extra terms to handle cosmic ray heating, the 
photoelectric effect, and radiation that is strongly frequency depen¬ 
dent. We replace equations 3 and 4 with 

4 (f) =-V F-V,:P-««.p{T.<-T.<)- 

acKip [T^ - Ta) , 


=-pV ■ V + acKgp - T^) - 

-A-gd Aline Ter “h Tpe + rH2,g , 

where equation 13 now describes the evolution of the specifle inter¬ 
nal energy of the gas, u, separately from the dust. The extra terms 
at the end of equation 13 are Agd which takes account of the ther¬ 
mal interaction due to collisions between the gas and the dust, Aune 
which is the cooling rate per unit volume due to atomic and molec¬ 
ular line emission. Ter which is the heating rate per unit volume due 
to cosmic rays, Tpe which is the heating rate per unit volume due to 
the photoelectric effect, and rH 2 ,g which is the heating rate per unit 
volume due to the formation of molecular hydrogen on dust grains. 
We also deflne the radiation constant a = 4aB /c and the dust tem¬ 
perature Td. We deflne the Rosseland mean gas opacity K,g (which 
is only important above the dust sublimation temperature) and the 
mean dust opacity na, which may be a Planck mean or Rosseland 
mean (see Section 3.3). 

We still need to determine the dust temperature. As done in 
most studies of the thermal structure of molecular clouds (e.g. 
Goldsmith 2001), we assume that the dust is in LTE with the 
total radiation held (i.e. that its temperature responds quickly to 
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any change). This replaces a dust thermal energy equation (which 
would also require the thermodynamic properties of the dust). 
Thus, we take 


pAisR + acHidP {Tr — Td) + Agd = 0, 


(14) 


where Aisr is the heating of the dust due to the interstellar ra¬ 
diation field (which is taken to be separate from the grey radia¬ 
tion field with temperature Tr). By rearranging equation 14 for Agd 
and eliminating this term from equation 13, it is easy to show that 
equations 12 and 13 reduce to equations 3 and 4 when Td = Tg 
and Aline ~ Agd ~ Ter ~ Tpe — rH 2 ,g ~ and where 
K,pp = (/^d + i^g)p and K,p = /^E. 

We note that the dust cooling rate 


Adust = acKdTd = 4crB QuBudu 
- j QjyBjydU 


CTB 




(15) 


= 47r 


where u is the frequency of the radiation, and Qu is the frequency- 
dependent dust absorption efficiency. Therefore when Agd = 0, 
and Tr = 0, and in the absence of extinction, the thermal balance 
with the interstellar radiation is obtained when 


Adust = Aisr = 47r 


/ 


^ tISRj 
dv ^ 


(16) 


where is the frequency-dependent ISR field. 

We note that equations 12 and 13 in Goldsmith (2001) appear 
to be incorrect. It is stated that the dust eooling rate is given by 




Adust — c / Ui/{Td) Ei{y^dv^ 


(17) 


where Uu{T) is the Planck energy density and where they took 

K,{iy) = 3.3 X cm“^, (18) 

with uq = 3.8x 10^^ Hz, and nH 2 is the number density of molec¬ 
ular hydrogen. However, we can write 

Adust — c J UdTd)K{v)di^ = 4Tr J BdTd)K{v)dv, 

3.3 X 10"^®nH247r 




I 


2hiy 


3.3 X 10"^®nH287r/!. 


/ 


q 2 ^hv/{kT^) — 

/_ I _ 

Qhu/{kT^) _ I 


dv, 


1 

di/, 


(19) 


3.3 X 10~^*^nH287r 


/< 


- 1 


ds, 


where h is Planck’s constant, k is Boltzmann’s constant, and we 
have written s = hu/{kTd). The integral can be performed over 
all frequencies, and can be evaluated as 


f 


^ds = r(6)C(6) = 5! C(6) = ^ ~ 122.08 (20) 


where r(n) = (n — 1)! is the Gamma function and the Riemann 
zeta function, C(2n) can be evaluated is 


C(2n) = 


(_l)»+iB2„(27r)" 

2(2n)! 


( 21 ) 


C(6) = 


(-l)‘‘B6(27r)® _ TT® 


2 ( 6 )! 


945 


( 22 ) 


where B 2 n is a Bernoulli number and Bq = 1/42. Thus, 


An.. ^122.08, 

= 4.22 X 10“^^nH2Td erg cm“^ s“ 


(23) 


This is a factor of 62 times larger (probably 2® = 64) than 
equation 13 in Goldsmith (2001) which states Adust = 6.8 x 
10“^^nH2Td. The effect of the greater dust cooling rate is to 
lower the dust temperatures by a factor of two for a given ISR 
field. We note that Glover & Clark (2012a) use Adust = 4.68 x 
10“^^nH2Td, which only differs from equation 23 by 10% (pre¬ 
sumably due to the choice of different dust opacities). 


2.4 Equations for heating and cooling terms in the ISM 

To calculate the extra heating and cooling terms appearing in equa¬ 
tions 12 to 14, we make the same approximations as those made in 
the past by many others who have studied the thermal structure of 
molecular clouds. 


2.4.1 Gas heating equations 

Following Goldsmith (2001), and Keto & Field (2005), we set the 
rate of energy transfer from cosmic rays into the gas as 

Fcr = 5 X 10~^®nH erg cm“^ s“^. (24) 


Note that Goldsmith (2001) and Keto & Field (2005) express their 
rates as functions of nH 2 = rip/2 for molecular gas. We will usu¬ 
ally refer to np throughout this paper, since when we allow for 
both atomic and molecular hydrogen the meaning of nH 2 may not 
be clear. 

The ISR is used to determine both the heating rate of the dust 
grains, and the photoelectric heating rate of the gas. In both cases, 
the ISR is attenuated due to dust extinction inside a molecular 
cloud. To describe the ISR, we use a slight modification to that de¬ 
scribed in detail by Zucconi, Walmsley & Galli (2001). They made 
approximate fits to the ISR given by Black (1994) using the sum of 
a power-law distribution and five modified blackbody distributions 
of the form 


tISR 



\ X ) ^ - 1 ’ 


(25) 


where A is the wavelength of the radiation, and the parameters Ap, 
Wi, and Ti are given in their Table B.l. Note that the mid-infrared 
term in their ISR is a power-law with a cut-off longward of 100 pm 
(which is mentioned in the main text of the paper, but not in their 
appendix). Although Zucconi et al. (2001)’s parameterisation does 
a good job of fitting Black’s ISR at most wavelengths, it does not 
provide a significant ultraviolet (UV) fiux which is necessary for 
photoelectric heating. To allow us to use one ISR parameterisation 
for both dust heating and the photoelectric effect, we add the ‘stan¬ 
dard’ UV background from equation 11 of Draine (1978) but only 
in the range hv — 5 — 13.6 eV. Note that to transform Draine’s 
parameterisation into the same units as equation 25 it must be mul¬ 
tiplied by h^v/eY. 

For the photoeleetric heating, we follow the preseription of 
Bakes & Tielens (1994), which has also been used by Young et al. 
(2004) and Keto & Caselli (2008), that 

Fpe = 1.33 X 10~^^e G(r)nH erg cm“^ (26) 


where, as described by Young et al. (2004), we have assumed that 


so that 
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the number density of nucleons Un ~ nn + 4 n(He) and that 
the gas is 25 percent helium by mass so that Un = 1.33 nn- The 
quantity G{r) is the ratio of the attenuated intensity of high energy 
(hu > 6 eV) photons to their unattenuated intensity 


G(r) = 


r47T POO 

/ / exp[—Tiy(r, cj)]di/dQ 

Jo JoeV 


p 4 : 7 T po 

Jo Joe 


(27) 




where Tu{r, u) is the frequency-dependent optical depth from r to 
the cloud surface along direction, u, and Q is solid angle. The effi¬ 
ciency factor e is a complicated function of the type of dust grains, 
radiation intensity, G{r), the temperature, and the electron number 
density. We use equation 43 of Bakes & Tielens (1994) 


_ 0.049 _ 

1 + 4 X 10-3 [G(r) Tg^^V(ne0PAH)]3'^3 

_ 0.037 (Tg/10^)°''^ _ 

1 + 2 X 10-‘^[G(r) Ty^/(ne0pAH)] 


(28) 


but with the adjustable parameter 0 pah that was introduced by 
Wolfire et al. (2003). We follow Wolfire et al. (2003) and use 
0PAH = 0.55 whereas the original equation of Bakes & Tielens 
(1994) had 0 pah = 1. For the conditions in cold molecular clouds 
the efficiency factor is nearly constant and can be approximated as 
e = 0.05 (note that in Keto & Caselli 2008, there is a typographi¬ 
cal error giving e = 0.5, and they also neglect the factor of 1.33 in 
Fpe). However, it is important to use the more complicated form at 
the lower densities of the warm and cold neutral mediums. Ordinar¬ 
ily the electron number density, rie, would come from a chemical 
model, but because we do not have such a model we need to pa- 
rameterise it. Using the results of Wolfire et al. (2003), in particular 
those displayed in their Fig. 10, we use the simple parameterisation 

Tie = nHmax(l x 10“^, min(l, O.OOS/nn)). (29) 


2.4.2 Gas cooling equations 


We develop a model for cooling in the diffuse ISM based on the 
results of the detailed model developed by Wolfire et al. (2003) 
(see also Glover & Mac Low 2007a). In the warm neutral medium 
(WNM), with a characteristic temperature of T ~ 8000 K, cool¬ 
ing is dominated by Lya emission from atomic hydrogen, electron 
recombination with small grains and polycyclic aromatic hydro¬ 
carbons (PAHs), and fine structure emission from atomic oxygen. 
At temperatures between those of the WNM and the cold neutral 
medium (CNM; T ^ 300 K), oxygen emission continues to con¬ 
tribute significant cooling but fine-structure emission from ionised 
carbon also becomes important and dominates in the colder parts 
of the CNM. Unlike Wolfire et al. (2003) and Glover & Mac Low 
(2007a), our aim is to achieve realistic temperatures without having 
to develop a detailed chemical model. Since we are not interested 
in regions of the ISM with temperatures greater than those of the 
WNM, we can produce a reasonable fit to the thermal behaviour 
by treating only the electron recombination, oxygen, and carbon 
emission. Following Wolfire et al. (2003) and Glover & Mac Low 
(2007a), we use the modified formula of Bakes & Tielens (1994) 

A„c = 4.65 X lO““0PAHT° ’^‘‘/nenH erg cm"® s“\ (30) 
where P = 0.74/Tg° °6®, and 


G{r)Tl/^ 

rie^PAH 


(31) 



Hh [cm-3] 

Figure 1. The equilibrium temperature as a function of the number density 
of hydrogen nuclei, nn, in the warm and cold neutral mediums achieved 
by considering cosmic ray and photoelectric heating (without extinction) 
balanced by cooling due to electron recombination, and fine structure emis¬ 
sion from carbon and oxygen (solid line; equations 24, 26, 30, 32, and 33). 
The dashed line gives the result obtained by Wolfire et al. (2003) for so¬ 
lar metallicity gas at a Galactic radius of 8.5 kpc. The dotted line gives the 
result from the model of Glover & MacLow (2007a). 


For the atomic oxygen fine-structure cooling, we use equation 
C3 of Wolfire et al. (2003) 

Aoi = 2.5xlO“^^nH exp ergcm“%"\ 

(32) 

Keto, Rawlings & Caselli (2014) calculate the cooling taking into 
account the abundance of atomic oxygen. 

For the atomic carbon fine-structure cooling, we use equation 
Cl of Wolfire et al. (2003) 

Aq+ = 3.15 X 10~^^nH exp erg cm(33) 

which assumes a carbon abundance relative to hydrogen of Ac — 
1.4 X 10“^. This is almost identical to equation of 2.67 of Tielens 
(2005), which was also used by Keto & Caselli (2008), except that 
in equation 33 the dependence on the degree of ionisation, Xi, has 
been neglected. Keto & Caselli (2008) assume that the degree of 
ionisation is equal to the abundance of C^ with respect to H 2 . 

At the low densities of the WNM and CNM, the equations for 
cosmic ray and photoelectric heating and cooling due to electron 
recombination and oxygen and carbon cooling can be used to cal¬ 
culate the equilibrium temperature as a function of density. This 
equilibrium curve resulting from our chosen parameterisations is 
shown in Fig. 1 as the solid line and compared to the results ob¬ 
tained by Wolfire et al. (2003) for solar metallicities at a Galactic 
radius of 8.5 kpc using a dashed line. Given that our model is much 
simpler than that of Wolfire et al. (2003) and that our main inter¬ 
est is in the star formation occurring in higher-density (nn ^ 100) 
molecular gas where other processes dominate, the level of agree¬ 
ment is satisfactory. Also plotted in the figure with a dotted line is 
the result obtained by Glover & Mac Low (2007a). We note that our 
model gives slightly lower temperatures than Wolfire et al. (2003) 
at densities nn ^ 20 cm“^, while the model of Glover & Mac Low 
(2007a) gives somewhat higher temperatures. 

For the molecular line cooling we follow Keto & Field (2005) 
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by using the parameterised cooling functions provided by Gold¬ 
smith (2001) for standard abundances given by 

Aline = a(Tg/10 K)^ erg cm“^ (34) 

where the parameters a and [3 are given in tables in Goldsmith 
(2001) as functions of nH 2 (we take nH 2 = nn/S). We allow 
the possibility of performing calculations both with standard abun¬ 
dances and with depleted abundances. In the former case, Table 2 
of Goldsmith (2001) gives cooling rates without allowing for de¬ 
pletion onto grains and covers the range nu 2 = 10^ — 10^ cm“^. 
We tried using linear interpolation of the logarithm of the cooling 
rates at the given points to compute the cooling rate at a particular 
value of nu 2 , but found that this tended to give discontinuities in 
the temperature gradients in molecular cloud cores. Instead we use 
three-point polynomial interpolation of the logarithm of the cooling 
rate (taking one point at the nearest tabulated value of nu 2 below 
the required density and two points above the required density, ex¬ 
cept for the highest densities where the reverse is done). This works 
well over the entire range of nH 2 except that it gives a ‘bulge’ in 
the cooling rate in the range nH 2 = 10^ — 10^ cm“^. To fix this, 
in this range, we calculate two three-point polynomial interpola¬ 
tions (one that has one point below the required value of nu 2 and 
two points above, and the other that has two tabulated values be¬ 
low and one above) and then use linear interpolation of the loga¬ 
rithm of these two values to obtain the cooling rate. At densities 
nH 2 < 10^, cm“^ we extrapolate the line cooling rate using the 
three-point polynomial derived from the first three points in the ta¬ 
ble. At densities nu 2 > 10® cm“®, we set the line cooling to zero 
(since it will be negligible compared to the other heating and cool¬ 
ing terms anyway). 

To allow for depletion, we use Table 4 of Goldsmith (2001) 
which gives cooling rates as functions of density over the range 
nH 2 = 10® —10® cm“® and as functions of a depletion factor rang¬ 
ing from 1 to 100. When depletion is allowed, for nH 2 <10^ cm“® 
we use the standard cooling rates divided by the depletion factor, 
and for nu 2 >10^ cm“® we use the standard cooling rates (since 
the cooling rates become less dependent on the depletion factor at 
high densities, according to Goldsmith 2001). In the intermediate 
density range we use bilinear interpolation in log-space of the val¬ 
ues in Table 4 that are given as functions of density and depletion. 
In the ranges nH 2 = 10^ — 10® cm“® and nH 2 = 10® — 10^ cm“® 
we further use linear interpolation in log-space to achieve smooth 
transitions between using the depleted cooling rates (from Table 4) 
and the depleted standard cooling rates (nH 2 = 10^ — 10® cm“®) 
or the standard cooling rates (nu 2 = 10® — 10^ cm“®). 


2.4.3 Dust heating and cooling equations 

For the heating rate of dust grains, we follow Zucconi et al. (2001) 
and Keto & Field (2005) and calculate the grain heating due to an 
incident radiation field whose intensity is attenuated by the optical 
depth averaged over 47r steradians upon grains of the absorption 
efficiency, Q^. Thus, 

pAn poo 

Aisr = / / QjyJl^^ex.p[-rjy{r,uj)]diydQ (35) 

Jo Jo 

where Tiy(r, u) is the frequency-dependent optical depth from r to 
the cloud surface along direction, cj, and we have deliberately omit¬ 
ted to normalise the integral over solid angle because of equation 
16. Zucconi et al. (2001) defines Qu in units of cm^ but in 
the above equation we redefine Qu to be in units of cm^ g“^ by 
dividing by /xmp, where rup is the proton mass. For simplicity we 


use the parameterisation of given by Zucconi et al. (2001). See 
section 3.3 for more information regarding the opacities we use. 

The thermal interaction between the gas and the dust depends 
sensitively on the distribution of grain sizes (Burke & Hollenbach 
1983). We use two different rates in Section 4, depending on which 
tests we are performing. For most of the calculations, we use the 
rate of Keto & Field (2005) 

Agd = 2.5 X 10“’^^ni[Tg^/2(Tg - Ta) erg cm“’^ s“\ (36) 


which is similar to that used by Goldsmith (2001) (our rate is a 
factor of v^/2 ^1.6 larger). However, this rate is more than an 
order of magnitude smaller than the rate given in equation 2.15 of 
Hollenbach & McKee (1989) 


Agd = 3.8 X 10"®®nHTg^^^(Tg - Td) 


X 1 


0.8 exp 



-3 -1 

erg cm s 


(37) 


assuming a minimum grain size of 0.01 fim. The last term in this 
equation takes account of the effects of the contributions of gas 
species other than protons and the effects of particle and grain 
charges. This is the rate used by Glover & Clark (2012a), so we 
use this rate for the tests in Section 4.4. 


2.4.4 Metallicity 

In all of the above equations, we have assumed standard abun¬ 
dances and gas-to-dust mass ratios. To allow the modelling of 
molecular gas with different metallicities, we assume that Qi, (and 
hence /^d), Agd, Aune, and Fpe all scale linearly with the metallic¬ 
ity. Thus, each of these quantities is multiplied by the factor Z/Z©. 
In doing so, we are also assuming that the grain properties are in¬ 
dependent of the metallicity. Note that Fcr does not depend on the 
metallicity. Also, the gas opacities, /^g, already explicitly include 
the metallicity dependence (see Section 3.3). 


2.5 Chemistry 

As mentioned in the introduction, our aim with this method is to 
develop a relatively simple model that captures the most important 
thermodynamic behaviour of the low-density interstellar medium, 
not to develop a full chemical model of molecular clouds. The 
equations in the previous section almost enable us to do this without 
modelling any chemistry at all. However, as we will see in Section 
4.4, we do need to have some model to treat the transformation of 
C+ into CO. At low-densities (uh ~ 10 — 10® cm“®) C+ is the ma¬ 
jor coolant, while at higher densities CO is the primary coolant until 
dust takes over (e.g. Glover & Clark 2012a). Furthermore, particu¬ 
larly at low metallicities, extra heating of the gas can be provided 
by the formation of molecular hydrogen so we need to treat the 
evolution of atomic and molecular hydrogen. 


2.5.1 Carbon chemistry 

Keto & Caselli (2008) provide a very simple model of the abun¬ 
dances of C^, neutral carbon, and CO, including the depletion 
of CO onto dust grains. Glover & Clark (2012a) show that the 
model significantly under-estimates the abundance of neutral car¬ 
bon, but this is only important in a narrow density range around 
nn ~ 10® cm“® and is unimportant for cooling (Glover & Clark 
2012b). Therefore, we implement the model of Keto & Caselli as 
stated in their paper, except as noted by Glover & Clark (2012a) 
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there is a typographical error in equation 5 of Keto & Caselli (2008) 
which should read 

CO ^ 6 X 10-i«nH2 

C+ 1.4x 10-iiGoexp(-3.2^v)’ ^ 

where Go is the ISR field in units of the Habing flux (Habing 1968) 
and we take Go = 1. 

The only significant extra quantity that we need to calculate 
in order to implement the chemical model of Keto & Caselli is the 
mean visual extinction 

(exp(-Av)) ^ ^ J exp(-Av)dQ, (39) 

where we take Ay = EQiy(V) where E is the column density in 
g cm“^ (so if the hydrogen is fully molecular the column density 
of H 2 is A^h 2 = E/[/imH]) and we take the frequency of visual 
light to be that of light with a wavelength of 550 nm. This is cal¬ 
culated at the same time as the integrals in equations 35 and 27 
(see Section 3.1.2). We use the resulting abundance to scale the 
fine-structure carbon cooling in equation 33, and we use the CO 
depletion factor that the model provides as the depletion value to 
use for the line cooling of Goldsmith (2001). 

2.5.2 Hydrogen chemistry 

For the hydrogen chemistry, we only consider H 2 formation on 
grains and dissociation due to cosmic rays and photodissociation. 
Omukai (2000) finds that grain formation dominates over gas phase 
formation (i.e. via H“ or HJ ions) for metallicities ^ 10“^. 

Glover (2003) compares gas-phase and grain-catalysed formation 
in detail and draws similar conclusions - for temperatures less than 
a few hundred Kelvin, H 2 formation on grains is expected to domi¬ 
nate for dust-to-gas ratios ^10 ^-10 ^of that found in the local 
ISM, over wide ranges of densities and levels of ionization. We 
neglect collisional dissociation of H 2 as at the temperatures we are 
dealing with in the low-density interstellar medium it is expected to 
be negligible compared to photodissociation and dissociation due 
to cosmic rays (c.f. the rates in Lepp & Shull 1983). 

In the above sections, when we have used nH 2 it has been in 
the context of fully molecular gas for which nH 2 = nn/S. How¬ 
ever, with the introduction of hydrogen chemistry, we need to dis¬ 
tinguish between atomic and molecular hydrogen. We therefore in¬ 
troduce the fraction of molecular hydrogen, xh 2 , which is equal to 
^H 2 /^H = 1/2 when the gas is fully molecular, and is equal to 
zero when the gas is fully atomic. 

For the formation rates of H 2 on grains, and the dissociation 
due to cosmic rays and photodissociation, we use the same model 
as Glover et al. (2010). The formation rate is given by equation 165 
in Table B1 of Glover et al. (2010), which is 

Rn2 = 3.0 X 10-^® T° Va/b {Z/Zq) cm® s”®, (40) 

where we have assumed that the rate scales linearly with the metal- 
licity due to the variation in the grain abundance, and where 

/a = [1.0 + 10^ exp(-600/Td)] , (41) 

/b = [1.0 + 0.04(Tg+Td)“® + 0.002Tg + 8xl0"®Tg]"\ (42) 
The number density of H 2 then evolves as 

= i^H2 [nu (1 — 2xh2)]^ . (43) 

form 

For the magnitude of the dissociation of molecular hydrogen 


dnH 2 

dt 


due to cosmic rays, we use the rate given in Table 2 of Bergin et al. 
(2004) 

= 1.2 X 10~^^xh 2 cm“^ s“^. (44) 

For photodissociation of molecular hydrogen we use the same 
model as Glover & Mac Low (2007a) and Glover et al. (2010), 
which is based on the work of Draine & Bertoldi (1996). We take 
the photodissociation rate of H 2 in optically thin gas to be 

= 5.6 X 10““ s“\ (45) 

where we have assumed a Draine (1978) UV interstellar radiation 
field. We take into account attenuation of the radiation due to dust 
absorption, and also self-shielding of the H 2 by line absorption due 
to other H 2 molecules. For the former we use 

/dust = (exp(-^v))®’’‘‘, (46) 


dnH 2 

dt 


which is simply a power of the mean visual extinction that is al¬ 
ready required for the carbon chemistry (equation 39). For the self¬ 
shielding we use 


/shield — 


0.965 0.035 


[- 8.5 X + , 

(47) 


where x — iVH 2/(5 x 10 ^^ cm“^), and 65 = 6 /( 10 ^ cm s“^), 
where b is the Doppler broadening parameter (which we take to be 
unity in the tests below). The magnitude of the fully shielded H 2 
photodissociation rate is then 


dnH 2 


dt 


pd 


-f^pd,o/dust/shield^H2 • 


(48) 


Once the total rate of change of molecular hydrogen has been 
obtained, the rate of change of the H 2 fraction is evolved as 

dxH2 _ J_ dnH2 .. Q. 

dt riH dt 


The formation of molecular hydrogen on dust grains releases 
approximately 4.5 eV of energy per molecule. A fraction of this 
will be radiatively lost while the remainder will heat the gas, with 
the relative fractions depending on the collisional de-excitation rate 
which depends on density. Following Glover & Mac Low (2007a), 
we assume that the heating rate of the gas is 


rH 2 ,g = 


7.2 X 10~^^ 

(1 + n cr /nu) 


dnH 2 

dt 


erg s ^ cm 

form 


where 


1 _ 1 — 2 xh 2 2 xh 2 

'blcr ricr,H ricr,H2 

logncr.H = 3.0 - 0.461 T 4 - 0.327 T|, 
log Tier,H2 = 4.845 — 1.3 T 4 + 1.62 T|, 


(50) 


(51) 


(52) 


with T 4 = Tg/10^ K and where the first equation, accounting for 
H 2 -H interactions, is an order of magnitude less than the value 
given by Lepp & Shull (1983), as recommended by Martin et al. 
(1996), and for H 2 -H 2 interactions the equation is taken from 
Shapiro & Kang (1987). 

In addition to heating due to molecular hydrogen formation, 
we investigated (using calculations similar to those presented in 
Section 4.4) the effect of including gas heating during H 2 pho¬ 
todissociation and due to UV pumping of H 2 (Glover & Mac Low 
2007a). However, we found that it is usually insignificant. The extra 



8 M. R. Bate & E. R. Keto 


heating only has a significant effect when the gas is initially highly 
molecular and has a low metallicity {Z % 0.1 Z©). Such circum¬ 
stances are highly unrealistic because the low extinction favours 
prior destruction of the H 2 . Even then, the heating only affects the 
outer, low-density (nn ^10^ cm“^) parts of the clouds (because 
of H 2 self-shielding) which are even less likely to be molecular 
than the inner parts of the clouds. Since we find this heating to be 
insignificant in realistic situations, and for the sake of simplicity, 
we do not include this effect. 


3 IMPLEMENTATION 

The calculations presented in this paper were performed using a 
three-dimensional smoothed particle hydrodynamics (SPH) code 
based on the original version of Benz (1990; Benz et al. 1990), 
but substantially modified as described in Bate, Bonnell & Price 
(1995), Whitehouse, Bate & Monaghan (2005), Whitehouse & Bate 
(2006), Price & Bate (2007), and parallelised using both OpenMP 
and the message passing interface (MPI). 

Gravitational forces between particles and a particle’s near¬ 
est neighbours are calculated using a binary tree. The cubic spline 
kernel is used and the smoothing lengths of particles are variable 
in time and space, set iteratively such that the smoothing length 
of each particle h — 1.2(rn /where m and p are the SPH 
particle’s mass and density, respectively (see Price & Monaghan 
2007, for further details). The SPH equations are integrated using 
a second-order Runge-Kutta-Fehlberg integrator (Fehlberg 1969) 
with individual time steps for each particle (Bate et al. 1995). To 
reduce numerical shear viscosity, we use the Morris & Monaghan 
(1997) artificial viscosity with varying between 0.1 and 1 while 
= 2av (see also Price & Monaghan 2005). 


3.1 Implementation of the new method 

In this section, we describe the detailed implementation of the 
method described in Section 2. The two main new elements are 
the implicit integration of equations 12 to 14, and the calculation 
of the dust attenuation which is required to compute the local ISR 
that appears in equations 27, 35, and the extinction in equation 39. 


з. 1.1 Solving the energy equations 

In solving the energy equations 12 to 14, we closely follow the im¬ 
plementation of the flux-limited diffusion method of Whitehouse 
et al. (2005) and Whitehouse & Bate (2006). The energy equations 
are very similar to those for the pure flux-limited diffusion method, 
except that the dust and gas now have different temperatures and 
there are some additional heating and cooling terms. Whitehouse 
et al. (2005) developed an implicit method for solving equations 3 
and 4 based on writing the equations in terms of the specific radi¬ 
ation energy, ^ = E/p, and the specific internal energy of the gas, 

и. For each SPH particle z, implicit expressions were derived for 

these two quantities at time step n + 1. These two equations were 
combined so as to solve for and The values from the 

previous time step, and were used as the initial guesses for 

and and Gauss-Seidel iteration was performed over all 
particles that were being evolved until the quantities converged to 
a given tolerance. The same basic method is used here. 

During each Gauss-Seidel iteration, the dust temperature must 
be solved for first. Equation 14 is solved directly for Td of each par¬ 
ticle based on its quantities from the previous time step or iteration. 


This is a root finding problem, for which Newton-Raphson iteration 
converges very quickly. It involves computing the lefthand side of 
equation 14 and its derivative with respect to Td, which is straight¬ 
forward except for the fact that is a function of Td. However, 
since is stored in a table as a function of Td (see Section 3.3), 
we simply use a numerical derivative to calculate d(/^d)/d(Td). 

Once the dust temperatures have been found, we solve for 
and in a similar manner to Whitehouse et al. (2005), 
but including the additional terms. However, we found it necessary 
to solve the equations in two different ways in the low-density and 
high-density regimes, refiecting the fact that equations 12-14 can 
be combined in two ways. In the low-density regime, when the gas 
and dust are poorly coupled, we solve 

'- 5 ©'-^ (53) 

+ pAisR + Agd, 

^ ) (54) 

Agd Aline T Ter T Tpe T rH2,g 5 


where we have obtained 53 by rearranging equation 14 to solve 
for the term acnap {T^ “ ^d) replaced this term in equation 
12. This works well in the poorly coupled regime because the Agd 
term is very small. Furthermore, in the low-density regime the gas 
temperature is usually Tg ^ 1000 K and hence K,g ^ 0. Thus, the 
two equations are almost uncoupled. However, in the well coupled 
regime, we find it better to solve 

- acKdp {T^ - Td) , 

p^ = -pV-v + acKgp {Tf - T© + acKAp {T* - T^) 

+ pAlSR — Aline + Ter + Tpe + rH2,g , 

(56) 


where we have obtained equation 56 by rearranging equation 14 to 
solve for Agd and replaced this term in equation 13. These equa¬ 
tions are almost identical to those solved by Whitehouse et al. 
(2005), except for the last four terms in equation 56, and these all 
tend to be very small when the gas and dust are well coupled. This 
form is better in the well-coupled regime because the Agd term in 
equations 53 and 54 can become very large when the gas and dust 
are well coupled even if the difference in the gas and dust temper¬ 
atures is very small (|Tg — Td| ^ 0.1 K). In this case, the Gauss- 
Seidel iterations can fail to converge and/or the radiation energy 
density can become negative. We use equations 53 and 54 when 
nii 2 {ZlZ@) < 10“ cm-3 and |Tr -Td| > 1 K, otherwise we use 
equations 55 and 56. 

We solve equations 53 and 54 or equations 55 and 56 implic¬ 
itly in a similar manner to Whitehouse et al. (2005). For equations 
55 and 56, we write 


= e?+di ^hc {piip^ - p,e; 
^ piPi 


.n+l\ '^Wjj 




— dtacua 


— dtacnd,i 




( 57 ) 


Pit. 




_ 7^4 

^ di 
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where 


b = 


Ai 


A,' 


Pi Pj 

for brevity, and 

—u^ + dt (pdV^) + dtac/^g,i 


(58) 




a 

V Cv,z ) _ 


- 10“^^dt 




~33^^^H2,z rpl /2 f 


pi 




P^ 
rcr..: 


10 K 


Tne.?; Z 


+ dt^ 

Pi Pi 


dt 


rH2,g,z 


(59) 


in which we have taken equation 36 rather than 37. The equations 
obtained when solving equations 53 and 54 are very similar (and 
simpler) to the above equations so we omit them for the sake of 
brevity. Note that for an ideal gas, the temperature of the gas is 
given by Tg = u/cv, where Cv is the specific heat capacity. In fact, 
our equation of state is considerably more complex (see Section 
3.2) and Cv is in fact simply the ratio of u/Tg which is obtained 
from a pre-calculated table of values that gives this ratio as a func¬ 
tion of u and p. Equation 59 includes the quantity pdV^ which is 
an explicit term (i.e. evaluated at time step n) for the work done 
on the gas and the thermal contribution from the artificial viscosity. 
In Whitehouse et al. (2005), these terms were included implicitly, 
but from Bate (2010) onwards we have used the explicit term here 
rather than the implicit term because we found empirically that us¬ 
ing the explicit term results in better energy conservation in star 
formation calculations. The explicit term for particle i is calculated 
as 


pdVi 


p. _^ 

Q 2 ^j'^ij ■ VzWij(/lz) 

'iri j 



rrij _ 

_ Ck,ij Vsig \ Vij 

Pij 


f-ifVij 


ViWifhi) ViWifhj) 

2 ni 20j 


(60) 


where the the second term is only applied between approaching 
particles (vij • rij < 0), the viscous signal velocity is Usig = 
Cs—2vij ‘Vij , and the average sound speed and densities of particles 
i and j are given by Cs = | (cs,z + Cs,j), and pij = \{pi + pj). re¬ 
spectively. The viscosity parameters evolve for each particle based 
on the method of Morris & Monaghan (1997), and their average is 

^ij — |(<^i + OLj). 

Within each Gauss-Seidel iteration, for each particle z, equa¬ 
tions 57 and 59 are solved simultaneously for and then 
This could be done exactly for the equations used by Whitehouse 
et al. (2005). Note, however, that equation 59 includes the quanti¬ 
ties Tg and (Tg,i)^^ in the collisional dust-gas term and the line 
cooling term, and the heating rate due to molecular hydrogen for¬ 
mation, rH 2 ,g, also involves the gas temperature. For a purely im¬ 
plicit solution of . the terms involving Tg,^ should be written 
as Tg i i. but this would introduce fractional powers of 

, making direct solution impractical. Instead, in equation 59 
we take Tg,i = u^/cv,z and rely on the fact that is updated in 
every Gauss-Seidel iteration so that if converges, so too does 
u2. Empirically, this seems to work well as long as large decreases 
in the gas temperature (which would result in a large decrease in the 


emission line cooling) are avoided from one iteration to the next — 
we limit ^ 0.8u-. 

Apart from these changes, the simultaneous solution for 
and then follows the method of Whitehouse et al. (2005) and 
will not be repeated here. Briefiy, in general, equations 57 and 59 
are combined to eliminate and produce a quartic equation for 
that can be solved analytically. When solving the implicit 
version of equation 54, we use Newton-Raphson iteration to obtain 
^n+i is obtained, is found simply by rearranging 

57. This is done for all active particles for each iteration, and the 
iterations are repeated until a desired tolerance is reached. 


3.1.2 Calculating the ISR attenuation 

It is necessary to calculate the local intensity of the ISR (attenu¬ 
ated due to dust extinction) and the mean visual extinction in order 
to evaluate the Aisr and Tpe terms and to model the chemistry 
(Section 2.5). This is potentially very difficult computationally be¬ 
cause it involves integrating over all angles at each point in the 
simulation. To make this computationally tractable, we use a sim¬ 
ilar method to the TREECOL method of Clark, Glover & Klessen 
(2012). They propose calculating the mass surrounding a particle 
in angular cones (and thereby calculating the column density and 
optical depths) using the same tree which is frequently used in SPH 
codes to calculate gravity and neighbouring particles. The algo¬ 
rithm essentially sums over the column densities of nearby par¬ 
ticles and more distant tree nodes that intersect the angular cone 
being considered. To cover the full 47r steradians uniformly, they 
used angular cones defined using the HEALPix^ library functions 
(Gorski et al. 2005). 

We implement a very similar method to that of Clark et al. 
(2012), with the following differences. The SPH code used by 
Clark et al. used an oct-tree. An oct-tree is constructed by a purely 
spatial decomposition of the particle distribution and the size and 
shape of each node in the tree is well-defined (in three dimensions it 
is a cube with a size that is purely a function of its level in the tree 
hierarchy). This allowed Clark et al. to roughly approximate the 
size and shape of each node projected onto the sphere as a square, 
thereby allowing them to easily calculate whether or not a node 
contributed to the mass lying in a particular angular cone. In turn, 
this allowed them to construct a method that conserved the total 
mass (i.e. average column density). However, the SPH code we use 
for this paper uses a binary tree based on a nearest neighbour de¬ 
composition (Benz et al. 1990). Thus, the nodes do not have well- 
defined shapes or sizes. We still implement a method that guaran¬ 
tees the average column density is exactly preserved, but for sim¬ 
plicity we assume that the projected shape of each node is a circle. 
Similarly, we assume that the projection of the solid angle covered 
by each HEALPix direction is a circle with area 47r/A^ steradians, 
where N is the chosen number of HEALPix directions. The binary 
tree defines a size for the nth node as 

f ^‘2 , mi , A 

Sn = max -^-ri2 + si,-^-ri2 + S2 , ( 61 ) 

\?ni+?n2 7711+7712 J 

where the subscripts 1 and 2 identify the two sub-nodes, rtj — 
\ri — rj|, and nrii and ri are the mass and position of the centre 
of mass of the sub-nodes, respectively. If the sub-node is an indi¬ 
vidual SPH particle, its size is zero in equation 61. As mentioned 
above, in general the tree nodes do not have well-defined shapes. 


1 


http://healpix.jpl.nasa.gov/ 
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For the purposes of calculating the column density we assume their 
projection is circular, but we are free to scale their effective radius 
by an arbitrary factor so that the radius of the circle is Vn = fsn, 
where / is expected to be of order unity. We explore the effects 
of choosing different values of / in Section 4.1 and conclude that 
/ = 0.5 is a good choice. If the node is an individual SPH parti¬ 
cle, we take the radius to be twice the particle’s smoothing length, 
En = 2/in. From a particular particle, p, at location, r^, the vec¬ 
tor between the particle and the node is Vnp = Tn — Vp. We then 
define the angular radius of the node as viewed from the particle 
by an = arctan(rn/'rnp). We denote the HEALPix direction as Z, 
which we treat as the centre of the HEALPix circle of angular ra¬ 
dius ttpix = y^4/A^ (the tt’s cancel). The angular distance between 
the centres of these two circles is cZ = arccos(rnp • Z/rnp). 

To calculate the contribution of each tree node or particle to 
the column density in each HEALPix area, we calculate the over¬ 
lapping area of the two circles. Eor two circles of radii r and R 
whose centres are separated by distance d this can be computed as 


A—r arccos 


+ R arccos 


+ r" - R^ 
2 dr 

d'^ + R^ - 


2dR 


— — \/(—d -|- r + -R)(d + r — IV)(^d — v R){d -|- r + IE ), 

(62) 


where we take r — apix and R — an. The contribution of the 
node to the average column density of the HEALPix area is then 
mn/{'Krn)A//N) where mnjij^r^ is the column density of 
the node and A/(47r/A^) is the fraction that is assigned to that 
HEALPix area. The computational efficiency can be improved by 
treating trivial cases separately. If Un Tupix ^ d there is no overlap. 
lfd + C^-pix ^ an then the HEALPix area is entirely covered by the 
node and the contribution is simply rUn/ (^rr^). If d an ^ apix 
then the node lies entirely within the HEALPix area and the contri¬ 
bution is mn/{4:7r/N). 

Once we have determined the average column density in a par¬ 
ticular HEALPix area, we can calculate the visual extinction Ay in 
that direction. We define Qv = Q^iV) at a visual wavelength of 
550 nm. Rather than performing the integrals appearing in equa¬ 
tions 27 and 35 during an SPH calculation, we pre-calculate tables 
of values. Eor equation 35, the attenuated dust heating in a single 
direction as a function of log(Av) 



Qviy) 
1.086 Qv 


di/, 


(63) 


where the factor of 1.086 = 2.5 log(e) appears because extinction 
A is related to optical depth r by A = —2.5 log(e“^). The integral 
is actually calculated in d(lnz/) rather than du. This table can be 
interpolated to find the heating rate for any particular value of the 
visual extinction Ay. The total heating rate at a location is then 
provided by averaging the contributions from all directions li as 


Aisr = ^ ^ X^^^{Avi), (64) 

i 


where the directions are chosen using the HEALPix library func¬ 
tions. The same is done to compute the quantity G(y) for the pho¬ 
toelectric heating (equation 27). We pre-calculate a table of values 


along a single line of sight containing 

9{M) = 


roo 

/ tISR 

/ exp 

JQeV 

Av Qu(y) 

dv 

1.086 Qv 

J 

roo 

1 

6eV 


( 65 ) 


The total value G(y) at a location (equation 27) is then provided by 
averaging the contributions from all directions U as 


1 ^ 

Gir) = 


( 66 ) 


Eor equations 39 and 46, we simply have 

1 ^ 

(exp(-^v)> = — y2exp(-^vi). (67) 

i 

To calculate the average column density of molecular hydrogen 
(A^h 2 ) which is required to evaluate the self-shielding of molec¬ 
ular hydrogen (equation 47) we use the same traversal of the tree 
that is used to calculate the extinction and G{r), but rather than 
calculate the full column density the contribution of each particle 
or tree node is multiplied by its molecular fraction xh 2 . 

It should be noted that even using the same tree that is used 
to calculate gravitational forces to estimate the extinction and other 
directional averages, obtaining the averages in many directions can 
still require a substantial amount of computation. Eurthermore, 
when the code is run in parallel using MPI with domain decom¬ 
position, the contributions of each domain to the column densities 
along each ray must be added together in order to calculate the 
quantities in each direction. This would not be necessary if one 
simply wanted the average column density to a particle (since the 
column density adds linearly), but because calculating the extinc¬ 
tion involves a nonlinear function, exp(—r), the column density in 
each direction must be computed separately. Thus, when running 
with MPI domain decomposition there is also a substantial mem¬ 
ory overhead which depends on the chosen number of directions. 


3.2 Equation of state and boundary conditions 

As in Bate (2009, 2012), the calculations presented in this paper 
were performed using radiation hydrodynamics with an ideal gas 
equation of state for the gas pressure p = pTglZ/fi. The thermal 
evolution takes into account the translational, rotational, and vibra¬ 
tional degrees of freedom of molecular hydrogen (assuming a 3:1 
mix of ortho- and para-hydrogen; see Boley et al. 2007). It also 
includes molecular hydrogen dissociation, and the ionisations of 
hydrogen and helium. The hydrogen and helium mass fractions are 
X — 0.70 and Y — 0.28, respectively. Eor this composition, the 
mean molecular weight of the gas is initially /x = 2.38. The contri¬ 
bution of metals to the equation of state is neglected. 

In previous calculations using the flux-limited diffusion 
method of Whitehouse et al. (2005), it was necessary to set mat¬ 
ter and radiation temperature boundary conditions. Eor calculations 
contained within a particular volume (e.g. the collapse of a molec¬ 
ular cloud core; Whitehouse & Bate 2006; Bate 2010, 2011; Bate 
et al. 2014) the matter and radiation temperatures were fixed on 
ghost particles outside the boundary at the initial (low) tempera¬ 
ture (e.g. 10 K). Eor calculations of clouds with free boundaries 
(e.g. Bate 2009, 2012, 2014), all particles with densities less than 
a particular value (typically 10“^^ g cm“^) had their matter and 
radiation temperatures set to the initial values (again 10 K). This 
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Figure 2. The dust temperature as a function of radius inside two uniform- 
density molecular cloud cores that are subject to an external interstellar ra¬ 
diation (ISR) field. The clouds have the same densities, but different masses: 
1 M 0 (upper distributions) and 10 M© (lower distributions). The exact so¬ 
lutions are plotted using the solid black lines. For the SPH calculations, 
a point is plotted for each SPH particle (2.6 x 10^ for each cloud). The 
different colours give the results obtained when the code uses 12 (red), 48 
(green), or 192 (blue) HEALPix directions to calculate the mean extinction. 
Dust within the clouds receives less of the external radiation due to extinc¬ 
tion and, therefore, is colder. The number of HEALPix directions used has 
little impact on the results, except in the very outer part of the clouds. The 
equilibrium temperature of dust that is subject to the full ISR is 17.2 K (hor¬ 
izontal dashed line), and the equilibrium temperature of dust that received 
exactly half of this radiation would be 15.3 K, which is lower by a factor of 
2^/® (horizontal dotted line). As expected, dust at the edges of the clouds 
is close to this temperature. 


matter was two orders of magnitude less dense than that in the ini¬ 
tial cloud and, thus, these boundary particles surrounded the region 
of interest in which the star cluster formed. In both cases, this essen¬ 
tially meant that the clouds were embedded in an external radiation 
field with this boundary temperature. 

The new method presented here eliminates the need for these 
arbitrary temperature boundary conditions. Now the temperature of 
the gas, dust and radiation are all set consistently at both high and 
low densities. 


3.3 Opacities and metallicity 

Whitehouse & Bate (2006), Bate (2009, 2010, 2011, 2012), and 
Bate et al. (2014) assumed solar metallicity gas and used the Rosse- 
land mean opacity tables of Pollack et al. (1985) for interstellar 
dust and, at higher temperatures when the dust is destroyed, the 
gas opacity tables of Alexander (1975) (the IVa King model) (see 
Whitehouse & Bate 2006, for further details). Bate (2014) per¬ 
formed calculations with varying opacities. He used the same dust 
opacities, scaled linearly in proportion to the metallicity, but re¬ 
placed the gas opacity tables with the metallicity-dependent tables 
of Ferguson et al. (2005) with X = 0.70. It is worth noting that 
one of the main conclusions of Bate (2014) was that the results of 



Figure 3. The dust temperature as a function of radius inside the 1 Mq 
uniform-density molecular cloud core that is subject to external ISR. The 
exact solution is plotted using the solid black line. The results from five 
SPH calculations are illustrated, and a point is plotted for each SPH parti¬ 
cle. The calculations each use 48 HEALPix directions to calculate the mean 
extinction, but different numbers of SPH particles are used to model the 
cloud - from top to bottom: 280 (black), 2608 (red), 2.6 x 10 ^ (green) 
2.6 X 10^ (blue), and 2.6 million (magenta). It can be seen that using fewer 
particles tends to result in warmer dust temperatures (i.e. an underestimate 
of the extinction), but the dust temperature appears to be slowly converging 
toward the exact solution as the number of particles is increased. As long as 
^3 X 10^ particles are used the maximum error in the mean temperature 
at any radius is ^ 1 K. The equilibrium temperature of dust that is subject 
to the full ISR is 17.2 K (horizontal dashed line), and the equilibrium tem¬ 
perature of dust that received exactly half of this radiation would be 15.3 K, 
which is lower by a factor of 2^/® (horizontal dotted line). As expected, 
with sufficient resolution, the dust at the edges of the clouds is close to this 
temperature. 


star formation calculations are very insensitive to the opacities that 
are used. 

In this paper, because the equations treating dust heating and 
photoelectric heating of the gas require integrals over the dust opac¬ 
ity as a function of frequency, we cannot continue simply to use the 
Rosseland mean opacities of Pollack et al. (1985). At low-densities 
(i.e. when the optical depth is low and the dust is essentially in 
thermal equilibrium with the interstellar radiation field), the Planck 
mean is more appropriate than the Rosseland mean and we desire 
consistency between the frequency-dependent dust opacity, Qu (z^), 
and the grey opacity, Kd appearing in equations 12 and 14. In other 
words, equation 15 should be satisfied. Therefore, we calculate 
Planck mean opacities directly from before the code be¬ 

gins to evolve the SPH calculation and the values are stored in a 
table as a function of dust temperature. Any frequency dependent 
opacities can be used in principle, but for simplicity, we use the 
parameterisation provided by Zucconi et al. (2001) of the opaci¬ 
ties from Ossenkopf & Henning (1994). We use these values for na 
whenever the dust temperature is less than 100 K. We assume that 
higher dust temperatures (Td > 100 K) will only be encountered 
at high densities (i.e. near protostars) when the gas and dust are 
optically thick and thermally well-coupled. In this case, Rosseland 
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Figure 4. The dust temperature as a function of radius inside the 1 Mq uniform-density molecular cloud core that is subject to external ISR. The exact solution 
is plotted using the solid black line. The results from ten SPH calculations are illustrated, and a point is plotted for each SPH particle. The calculations each use 
48 HEALPix directions to calculate the mean extinction and 2.6 x 10^ SPH particles, but the tree-opening parameter and the effective sizes of the tree nodes 
are varied. In the left panel, the tree-opening parameter is 0.5, while in the right panel it is 0.25. In each panel the effective size of the tree nodes decreases from 
top to bottom, using scaling factors of / = 2.0 (red), 1.0 (green), 0.5 (blue), 0.1 (magenta), and 0.01 (black). It can be seen that using very large tree-opening 
angles and/or effective node sizes results in warmer dust temperatures (i.e. an underestimate of the extinction). If the effective node size is too small the scatter 
increases. With a tree-opening parameter of 0.5, an effective node size of half the actual node size (i.e. / = 0.5) gives the best trade off between accurately 
computing the mean radial temperature distribution and minimising the scatter. The equilibrium temperature of dust that is subject to the full ISR is 17.2 K 
(horizontal dashed line), and the equilibrium temperature of dust that received exactly half of this radiation would be 15.3 K, which is lower by a factor of 
2^/® (horizontal dotted line). As expected the dust right at the edges of the clouds is close to this temperature. 


mean opacities are appropriate and we use those of Pollack et al. 
(1985) as we have in our earlier calculations mentioned above. 

For the gas, we continue to use the Rosseland mean opacities 
of Ferguson et al. (2005) for since the gas opacity only becomes 
important in the highly optically-thick regions inside, or very near 
to, a star, for which the Rosseland mean is appropriate. 

4 CALCULATIONS 

4.1 Testing the method for calculating ISR attenuation 

One of the tests that Clark et al. (2012) used for the TREECOL 
method was to calculate the temperature structure of two uniform- 
density spherical clouds when subjected to a Black (1994) ISR 
field. They took densities of 10“^^ g cm“^ for clouds with masses 
of 1 M© and 10 M© and used 2.6 x 10^ SPH particles. We per¬ 
form the same tests here. Eor this test, we set the internal radiation 
to zero (i.e. = Tr = 0) and we turn off the coupling between 

the gas and the dust (i.e. Agd = 0) so that equation 14 is only 
solving for the equilibrium temperature of the dust subject to the 
attenuated external ISR. The ISR is as discussed in Section 2.4.1, 
except that we exclude the additional UV fiux which is important 
for photoelectric heating. We have calculated the exact radial tem¬ 
perature distribution by analytically calculating the column density 
along 192 HEALPix directions (using 48 provides an almost iden¬ 
tical result) and used these values to calculate the attenuated dust 
heating (equations 63 and 64) and the equilibrium dust temperature 
as functions of radius. 

In Pig. 2 we give the distribution of the equilibrium dust tem¬ 
perature as a function of radius for each of the two clouds. One 


point is plotted for each SPH particle. Different colours give the 
results obtained when the code uses 12, 48, or 192 HEALPix di¬ 
rections to calculate the mean extinction. The number of HEALPix 
directions used has little impact on the results, except in the very 
outer part of the clouds when using only 12 directions. However, 
using more directions takes longer to compute, so from this point 
on we use 48 HEALPix directions. The equilibrium temperature of 
dust that is subject to the full ISR is 17.2 K for our chosen ISR. By 
considering equation 23, we expect a dust particle that is subject to 
one hemisphere of heating to have an equilibrium temperature that 
is a factor of 2^^® lower (i.e. 15.3 K) and, as expected, this is es¬ 
sentially equal to the temperature at the edge of each of the clouds. 
However, within the clouds we find that the dust temperatures are 
up to IK warmer than given by the exact solutions (solid black 
lines). We note that the central temperatures in our clouds are in 
good agreement with those found by Clark et al. (2012), but they 
obtained substantially lower temperatures at the outer edges of their 
clouds (13 — 14 K). The reason for this is not clear since, as just 
mentioned, our edge temperatures are what we would expect. Their 
ISR may have been different to ours, although both are nominally 
based on that of Black (1994). 

In Pig. 3, we plot the dust temperatures from five calculations 
of the 1-M© cloud that each use 48 HEALPix directions, but that 
have different numbers of SPH particles: 280, 2608, 2.6 x 10^, 
2.6 X 10^, and 2.6 million. It can be seen that the dust temperature 
tends to be over-estimated when a smaller number of particles is 
used, but the dust temperature appears to be slowly converging to¬ 
ward the exact solution as the number of particles is increased. The 
error in the dust temperature is approximately halved each time the 
number of SPH particles is increased by an order of magnitude (i.e. 
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the error decreases at about the same rate that the linear spatial res¬ 
olution increases). Comparing the 2.6 x 10^ particle calculation 
with the exact solution, the maximum difference between the mean 
temperatures at any radius is only 1 K. Clark et al. (2012) did not 
discuss how their method behaves with different numbers of parti¬ 
cles. 

We used a cubic lattice SPH particle distribution to generate 
the results presented in Figs. 2 and 3. We have also tried a random 
particle distribution. Using 2.6 x 10^ randomly-placed SPH parti¬ 
cles to model the l-M© cloud we obtained an almost identical mean 
radial temperature distribution to the case that used a cubic lattice, 
but with a temperature scatter of up to ±0.2 K. However, given that 
this error is similar to the systematic difference between the tem¬ 
perature distributions obtained from the 2.6 x 10^ and 2.6 million 
particle cubic lattice calculations, we conclude that the accuracy of 
the solution depends more on the number of particles used than on 
the details of how they are distributed. 

In Fig. 4, we investigate the effects on the calculation of the 
dust temperature of varying the tree-opening criterion and the ef¬ 
fective radius of the tree nodes, the latter of which is parameterised 
by the factor / (see Section 3.1.2). During the walk through the 
tree to calculate gravity, nodes are opened if Sn/^np is larger than 
a critical value, usually taken to be 0.5. We perform calculations 
using a tree-opening parameter of 0.5 and 0.25 (the latter of which 
means more nodes are opened during the tree walk). For each case, 
we perform five SPH calculations of the l-M© cloud that each use 
48 HEALPix directions and 2.6 x 10^ SPH particles, but we vary 
the effective radius of the node that is used when calculating the 
column density. We take factors of / = 2.0, 1.0, 0.5, 0.1, and 0.01. 
It can be seen from Fig. 4 that using very large tree-opening an¬ 
gles and/or effective node radii results in warmer dust temperatures 
(i.e. an underestimate of the extinction). We also note that almost 
identical results are obtained when the product of the opening pa¬ 
rameter and the factor / is a constant (e.g. using a tree-opening 
parameter of 0.5 and / = 1 gives almost identical results to using 
a tree-opening parameter of 0.25 and / = 2). The calculation of 
the extinction becomes very poor when the product of these two 
quantities exceeds 0.5. If the effective node size is too small, the 
calculations provide reasonable mean temperatures at a given ra¬ 
dius, but the scatter increases (e.g. the black points in the left panel 
that were produced using / = 0.01). With the typical tree-opening 
parameter of 0.5, setting the factor / = 0.5 gives the best trade off 
between accurately computing the mean radial temperature distri¬ 
bution and minimising the scatter. 

The typical resolutions employed in SPH calculations of star 
cluster formation vary from ~ 10^ particles per solar mass (Bon- 
nell. Bate & Vine 2003; Bonnell et al. 2011) to ~ 10^ particles per 
solar mass (Bate, Bonnell & Bromm 2003; Bate 2012). The results 
of this test suggest that the typical errors in the dust temperature 
due to finite numbers of SPH particles and HEALPix directions 
and variations of the tree parameters should be ^ 1 K as long as 
^3x10^ particles per solar mass are used and large effective node 
radii are avoided. 

4.2 Equilibrium dust and gas temperatures 

More realistic than using a uniform-density sphere is to use Bonner- 
Ebert spheres. To allow comparison with earlier work, we use the 
same 5-M0 Bonner-Ebert spheres that were used by Keto & Eield 
(2005) — a (marginally) subcritical case with a central density of 
nH 2 = 10^ cm“^, and a supercritical case with a central density 
of nH 2 = 10® cm“^. The former case has a ratio of the inner to 



Figure 5. Plots of the SPH molecular hydrogen number density, nii 2 , as 
functions of radius inside the subcritical (blue) and supercritical (red) 5 M© 
Bonner-Ebert spheres. A point is plotted for each SPH particle. The ‘spikes’ 
in the density for the supercritical case are an artefact of the particles being 
set up on a radially-deformed cubic lattice. 

outer density of 14, while the ratio of the latter is 3000. Their radii 
are 0.223 and 0.265 pc, respectively. In these tests, we assume the 
hydrogen is fully molecular. We model both clouds with 3 x 10® 
SPH particles. 

The SPH particles were set up using a cubic lattice which 
was then deformed radially to obtain the required cumulative radial 
mass profile. In Pig. 5 we provide the SPH density as a function of 
radius for of each SPH particle for both cores. The small ‘spikes’ in 
the density in the supercritical case are due to the lattice structure 
affecting the density calculation of some particles in this case with 
a very strong density gradient. They do not appear in the subcriti¬ 
cal case, which has a much shallower density profile. In both cases, 
there is also some ‘noise’ in the density near the outer boundary 
(which is made of refiected ghost particles). 

In this section, we investigate the effects of each of the phys¬ 
ical heating and cooling mechanisms on the temperatures of the 
gas and dust. We begin by including only cosmic ray heating and 
molecular line cooling for the gas, and the dust is taken to be in 
thermal equilibrium with the ISR. We then add in the effects of 
gas-dust collisions, photoelectric heating of the gas, and cooling 
from C^. Por the first two cases, we also investigate the effects of 
including and excluding the UV flux. 

4.2.1 Cosmic ray heating, molecular line cooling, and dust 
radiative equilibrium 

In Pig. 6, we plot the dust and gas temperatures as functions of 
radius inside the two Bonner-Ebert spheres. These calculations in¬ 
clude only cosmic ray heating and undepleted molecular line cool¬ 
ing of the gas, and the dust is in thermal equilibrium with the ex¬ 
ternal ISR. There is no gas-dust collisional coupling, photoelectric 
effect, or cooling due to C^, oxygen, or recombination lines. There 
are two calculations for each case — one that excludes the UV con¬ 
tribution to the ISR, and one that includes it. 

In the outer parts of the cores the dust temperature exceeds the 
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Figure 6. The dust (upper) and gas (lower) temperatures as functions of radius inside the subcritical (left) and supercritical (right) 5 Mq Bonner-Ebert spheres. 
The ISM physics includes cosmic ray heating and molecular line cooling of the gas, and the dust is in thermal equilibrium with an external ISR that excludes 
(green) or includes (red) the UV contribution (see Section 2.4.1). The dust temperature drops as the extinction increases at smaller radii. The gas temperature 
rises because the line cooling becomes less effective at higher densities. The UV ISR only affects the dust temperature in the low-extinction outer parts of the 
cloud. The small ‘spikes’ in the gas temperature in the right plot at radii r 0.03 — 0.16 pc are due to the artefacts in the calculation of the SPH density due 
to the particles being set up on a radially-perturbed cubic lattice. 




Radius [pc] Radius [pc] 

Figure 7. The dust (upper) and gas (lower) temperatures as functions of radius inside the subcritical (left) and supercritical (right) 5 Mq Bonner-Ebert spheres. 
Red and green points (mostly obscured) are the same as in Fig. 6. For the blue (without the UV ISR) and magenta (with the UV ISR) points, the ISM physics 
is exactly the same, except that it includes thermal collisional coupling between the gas and the dust. The black solid lines give the results obtained by Keto & 
Field (2005) without including the UV ISR. The effect is that at high densities, the gas essentially adopts the dust temperature. The small ‘spikes’ in the gas 
temperature in the right plot at radii r ^ 0.03 — 0.16 pc are due to the artefacts in the calculation of the SPH density due to the particles being set up on a 
radially-perturbed cubic lattice. 
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Figure 8. The effect of photoelectric heating and C+ cooling of the gas on the dust and gas temperatures as functions of radius inside the subcritical (left) 
and supercritical (right) 5 Mq Bonner-Ebert spheres. The dust temperatures are essentially unaffected by the extra gas heating and cooling processes. The 
blue points are the same as in Fig. 7, where the gas is subject to cosmic ray heating, molecular line emission, and collisions with the dust, while the dust is 
subject to heating from the ISR (including the UV), thermal emission, and collisions with the gas. The green points include the same physical effects as the 
blue points, but with the addition of photoelectric gas heating (equation 26). The photoelectric heating has an enormous effect on the gas temperature in the 
outer, low-density parts of the clouds. However, when the C+ cooling is included (red points), the effect of the photoelectric heating on the gas temperature is 
greatly reduced. 




Radius [pc] Radius [pc] 

Figure 9. The effect of CO depletion on the gas on the dust and gas temperatures as functions of radius inside the subcritical (left) and supercritical (right) 
5 Mq Bonner-Ebert spheres. The red points are the same as in Fig. 8, where the gas is subject to cosmic ray and photoelectric heating, molecular and C+ 
line emission, and collisions with the dust, while the dust is subject to heating from the ISR (including the UV), thermal emission, and collisions with the 
gas. The black points include the same physical processes as the red points, but also allow for CO depletion. This reduces the molecular line cooling at high 
densities which raises the gas temperature slightly at densities nH 2 = 10^ — 10^ cm“^. At higher densities the depletion has little effect because the cooling 
is dominated by dust emission through collisions with the dust. 
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Figure 10. The effects of changing the strengths of the cosmic ray (left) and photoelectric (right) heating of the gas on the dust and gas temperatures as 
functions of radius inside the supercritical 5 Mq Bonner-Ebert sphere. In both panels the red points are the same as in the black points in the right panel 
of Fig. 9. The other points include the same physical effects as the red points, but with one order of magnitude less (green) or more (blue) cosmic ray (left) 
or photoelectric (right) heating. Changing the level of cosmic ray heating has a large effect on the gas temperature at intermediate densities, but not at high 
densities (where the gas and dust are well coupled) or low densities (where photoelectric heating dominates). Changing the levels of photoelectric heating has 
a large effect on the gas temperature only in the outer regions of the core. The dust temperatures are unaffected in either case. 


gas temperature as the local ISR is only weakly attenuated by the 
dust, and the cosmic ray flux and line cooling are such that the low- 
density equilibrium gas temperature is ^ 10 K. Including the UV 
contribution to the ISR boosts the dust temperature by up to 2 K in 
the outer parts, but the total energy in the UV flux is small com¬ 
pared to the total ISR flux and the UV does not penetrate very far 
into the cores. In the inner parts of the cores, the gas temperature 
rises as the effectiveness of the line cooling decreases (Goldsmith 
2001). On the other hand, the dust temperature decreases, particu¬ 
larly in the supercritical case, because the dust extinction attenuates 
the ISR. 


4.2.2 The effects of gas-dust collisions 

In Fig. 7, we include the same physical processes as in the previ¬ 
ous section, but we now also turn on the transfer of thermal energy 
between the gas and the dust due to collisions. The strength of this 
interaction depends on the square of the density (equation 36), so 
we expect it to have little impact at low densities, but a signiflcant 
impact as the density increases. Indeed, this can be clearly seen in 
Fig. 7, in which we plot the same calculations as in Fig. 6 for ref¬ 
erence, but we now also include the case with gas-dust coupling in 
blue (without the UV ISR) and magneta (with the UV ISR). The 
gas-dust thermal coupling has no signiflcant effect on the tempera¬ 
tures in the subcritical core because the densities are too low. In the 
supercritical case, the dust temperature is unaffected, but the gas 
temperature at r < 0.04 pc (densities nH 2 ^ 2 x 10^) is ‘dragged 
down’ by the interaction with the cold dust so that at the centre the 
gas and dust temperatures are both 7.5 K. 

Our numerical results are in reasonable agreement with those 
presented in Figs. 4 and 5 of Keto & Field (2005), who use the same 
models for the ISM physics as we use here, without the UV ISR 


(black solid lines in Fig. 7). The main difference is that the gas tem¬ 
peratures at low densities are slightly lower in Keto & Field (2005). 
The gas temperature at low density is set by the balance between 
cosmic ray heating and the line cooling rates. The latter is depen¬ 
dent on the method used to interpolate from the tables of Goldsmith 
(2001). Since the cosmic ray heating and line cooling rates are both 
simple functions of nu (equations 24 and 34), it is easy to solve for 
the central gas temperature of 12.4 K in the subcritical case (where 
the gas-dust coupling is negligible). The slightly lower gas temper¬ 
ature 11.5 K) obtained by Keto & Field (2005) was found to be 
due to less accurate interpolation used by Keto & Field of the line 
cooling functions of Goldsmith (2001). Thus, the different interpo¬ 
lation of the line cooling produces the different gas temperatures. 


4.2.3 The effects of photoelectric heating and carbon chemistry 

In Fig. 8, in addition to the other physical processes, we turn on 
photoelectric heating of the gas due to UV ISR photons liberat¬ 
ing electrons from dust grains (equation 26). The resulting gas and 
dust temperatures (green) are compared to those without the pho¬ 
toelectric heating (blue). The dust temperature is essentially un¬ 
changed, but the gas temperature in the outer parts of the cores, 
at densities such that the gas and dust are thermally decoupled 
(nH 2 ^10^ cm“^), is much hotter than without photoelectric heat¬ 
ing. In fact the photoelectric heating is so strong that the cosmic ray 
heating is irrelevant in the low-density gas. 

Whereas Keto & Field (2005) did not consider photoelectric 
heating, as discussed in Section 2.5, Keto & Caselli (2008) included 
photoelectric heating and a simple carbon chemistry model which 
allowed them to include both cooling from at low densities and 
treat the depletion of CO at high densities. Therefore, in Fig. 8, 
we keep the photoelectric heating on, but also add cooling from 
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Figure 11. The effects of changing the strength of the interstellar radiation 
(ISR) field on the dust and gas temperatures as functions of radius inside 
the supercritical 5 Mq Bonner-Ebert sphere. The red points are the same 
as the black points in the right panel of Fig. 9. The other points include the 
same physical effects as the red points, but with one order of magnitude 
less (green) or more (blue) ISR heating. Changing the level of ISR by an 
order of magnitude has a 50% effect (10^/®) on the dust temperatures. At 
high densities the different dust temperature affects the gas temperature in 
a similar manner, but at low-densities the gas temperature is unaffected as 
it is set by the photoelectric heating and C+ cooling. 

(red). It can be seen that the cooling dramatically lowers 
the temperatures in the outer parts of the cloud; the temperature is 
still somewhat higher than it was without either the photoelectric 
heating or the cooling, but much less than if cooling is 
omitted. We note that recombination and oxygen cooling are both 
insignificant at these densities and that the temperatures are essen¬ 
tially identical whether they are included or not. 

Finally, in Fig. 9 we turn on the effects of CO depletion based 
on the equilibrium prescription of Keto & Caselli (2008). This re¬ 
sults in only slightly higher gas temperatures at intermediate densi¬ 
ties (nH 2 = 10^ — 10^ cm“^) where the densities are still too low 
for dust cooling to dominate, but the densities are high enough that 

is ineffective. 

4.2.4 The effects of changing the ISM parameters 

The above sections used the diffuse ISM model parameters de¬ 
scribed in Sections 2.3 and 2.4. In this section, we investigate the 
dependence of the gas and dust temperatures on the strengths of 
the various external heating sources. We only use the supercritical 
Bonner-Ebert sphere in these tests because this has the larger range 
of densities. In the left panel of Fig. 10 we explore the effects of 
changing the cosmic ray heating rate (equation 24) by an order of 
magnitude in each direction. In the right panel of Fig. 10 we explore 
the effects of change the gas photoelectric heating rate (equation 
26) by an order of magnitude in each direction. It can be seen that 
over particular density ranges, either of these changes the gas tem¬ 
perature by approximately half an order of magnitude in each direc¬ 


Figure 12. The effects of changing the metallicity on the dust and gas tem¬ 
peratures as functions of radius inside the supercritical 5 M© Bonner-Ebert 
sphere. Photoelectric heating has been included. From top to bottom, the 
calculations have metallicities of 1/100 (magenta points, top curves), 1/10 
(blue points), solar metallicity (red points), and 3 times solar metallicity 
(green points). The red points are the same as in the black points in the right 
panel of Fig. 9. The dust temperatures are greater with lower metallicities 
because the ISR is less attenuated by extinction. The gas temperatures are 
greater with lower metallicities because the UV ISR is less attenuated by 
extinction (leading to stronger photoelectric heating) and also because the 
gas emission line cooling is reduced. 


tion because the gas line cooling scales roughly as the square of the 
gas temperature for the molecular lines (Goldsmith 2001). The gas 
temperatures in the outer parts of the cloud are much more strongly 
affected by photoelectric heating than cosmic ray heating. On the 
other hand, at number densities nH 2 ^10^ cm“^ strong cosmic 
ray heating can raise the gas temperature significantly above the 
dust temperature whereas the photoelectric heating has no effect 
this deep within the core. 

In Fig. 11, we probe the effects of changing the ISR field 
(equation 25) by an order of magnitude in each direction, excluding 
the UV contribution which is held fixed. The calculations include 
photoelectric heating of the gas and cooling. The increased ISR 
primarily affects the dust temperature, but because the dust emis¬ 
sion scales as the sixth power of its temperature (equation 23), 
this only has a 50% effect on the dust temperatures. Deep within 
the core where the gas becomes thermally coupled to the dust, the 
warmer or cooler dust also results in warmer or cooler gas, respec¬ 
tively. 

In Fig. 12, we investigate the effects of changing the metallic¬ 
ity of the core, performing additional calculations with metallici¬ 
ties of 1%, 10%, and 3 times solar. It can be seen that the highest 
metallicity results in the coldest temperatures for both the gas and 
the dust due to the combination of several effects. The increased 
metallicity increases the ISR attenuation, decreasing the dust tem¬ 
perature inside the core. Similarly, the UV is prevented from pen¬ 
etrating as far into the core, reducing the effects of photoelectric 
heating of the gas, and the gas line emission is enhanced. Finally, 
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Figure 13. The evolution of the central gas and radiation temperatures as 
functions of the maximum (central) density during a radiation hydrodynam- 
ical calculation of the collapse of an unstable spherical 5 Mq Bonner-Ebert 
sphere. The radiation temperature is always less than or equal to the gas 
temperature in this graph. After the central density exceeds 10“^^ g cm“^, 
the central dust temperature is indistinguishable from that of the gas. The 
black solid lines give the results from a calculation using the method de¬ 
scribed in this paper, while the red dashed lines and blue long-dashed 
lines give the results using the pure flux-limited diffusion (FED) method of 
Whitehouse & Bate (2006). The diffuse ISM model used in the former cal¬ 
culation determines the initial gas and dust temperatures, which vary with 
radius, and the internal radiation temperature is arbitrarily set to a low initial 
value (1 K). In the FED calculations the initial temperatures of the matter 
(gas and dust) and internal radiation are arbitrary and we set them to 14 K 
(red dashed) and 10 K (blue long-dashed). As expected, it is mainly the 
low-density evolution that is affected by including the diffuse ISM model. 

there is an increase in the effectiveness of collisional thermal cou¬ 
pling between the gas and the dust. With the lowest metallicity, the 
dust almost becomes isothermal since there is almost no attenua¬ 
tion of the ISR, and the gas is warmer because the photoelectric 
heating is strong while emission line cooling is weak. 

4.3 The evolution of a collapsing Bonner-Ebert sphere 

In this section, as opposed to the static tests from the previous sec¬ 
tions, we report the results from a hydrodynamical calculation of 
the collapse of a Bonner-Ebert sphere and compare them to the 
results obtained using the pure flux-limited diffusion method of 
Whitehouse et al. (2005) and Whitehouse & Bate (2006). We be¬ 
gin with a somewhat unstable Bonner-Ebert sphere, similar to those 
studied in the previous section. We choose a S-M© core with an in¬ 
ner to outer density contrast of 20 and a radius of 0.1 pc. We model 
the core with 3 x 10^ SPH particles and use spherical reflective 
boundary conditions (modelled using ghost particles). 

In Eig. 13 we plot the central temperatures as functions of the 
central (maximum) density for the two methods. All of the initial 
temperatures in the pure flux-limited diffusion calculations are ar¬ 
bitrary. As is standard practice in such calculations, we set the mat¬ 
ter and radiation temperatures to be in thermal equilibrium initially. 
We perform two separate calculations and with two different initial 
temperatures of 10 K and 14 K. For the calculation performed using 
the method described in this paper, the initial gas and dust temper¬ 
atures are not arbitrary — they are set by our model of the diffuse 
ISM. The initial dust temperature increases from 8 K in the centre 


to 17 K at the outer edge of the molecular cloud core. The initial gas 
temperature is slightly warmer than the dust at the centre (10 K) and 
cooler at the outside 15 K). The initial central dust and gas tem¬ 
peratures are quite similar because the central density of the core is 
quite high (p = 6 x 10“^^ g cm“^ or nH 2 = 1.5 x 10^ cm“^) 
so that there is signiflcant thermal coupling between the dust and 
the gas (we use equation 36 in these calculations). As the collapse 
proceeds, the central dust and gas temperatures quickly converge, 
so we do not plot the central dust temperature in Fig. 13. The in¬ 
creases in extinction and density during the collapse mean that both 
the dust and gas temperatures drop to a minimum of 6.5 K when the 
central density is 10“^® g cm“^ before the compressional heating 
starts to increase the central temperature again. The initial radia¬ 
tion temperature should be set to a low value, because the heating 
of the initial conditions is via external radiation rather than inter¬ 
nal radiation. As long as a small value is chosen the exact value is 
unimportant, but we wish to avoid setting it to zero so as to avoid 
numerical problems. We chose a value of 1 K. 

As may be expected. Fig. 13 shows that it is mainly the low- 
density phase of evolution of the clouds that differs between the two 
methods (densities ^ 10“^^ g cm“^). Once the central regions of 
the core become optically-thick to infrared radiation (the so-called 
opacity limit for fragmentation; Eow & Eynden-Bell 1976; Rees 
1976) all of the temperatures converge. The evolution obtained 
with the pure flux-limited diffusion calculation with the lower ini¬ 
tial temperature (10 K) is closest to the evolution obtained using 
the method presented here, presumably because the central tem¬ 
perature in this model is closest to that given by the ISM initial 
conditions. 

In Fig. 14 we plot temperature profiles as functions of radius 
(left panel) and density (right panel) at various points during the 
collapse. Gas is shown with solid lines, dust with dotted lines, and 
radiation with dashed lines. At either a given radius or density the 
temperatures tend to rise during the collapse as the protostar emits 
more radiation. The formation of the stellar core at late times (ini¬ 
tial radius 0.01 AEl and temperature > 10^ K) is clearly visible 
in the left panel. Except in the outer parts of the core, these profiles 
are very similar to those obtained using pure flux-limited diffusion. 
The main difference in the outer parts of the core (radii ^10^ AU) 
is that the gas, dust, and radiation all have different temperatures 
with the new method, whereas with pure flux-limited diffusion all 
are very close to the initial arbitrary temperature (e.g. 10 K; see 
Whitehouse & Bate 2006). 

After stellar core formation as radiation works its way out¬ 
wards from the centre the gas, dust, and radiation temperatures 
become inverted at intermediate radii in the core (radii 400 — 
3000 AEl; Fig. 15). Before this point the gas temperature is the 
hottest and the radiation temperature the coldest at these radii. 
However, the radiation emitted by the proto star quickly dominates 
that from the external radiation held at these radii. The dust adopts a 
new thermodynamic equilibrium with the combined radiation held, 
being heated to temperatures of 15 — 40 K, but the gas which is 
poorly coupled to either the radiation or the dust at these low den¬ 
sities remains at temperatures of 8 — 12 K. Such effects are impos¬ 
sible to capture unless the thermal evolution of the gas and dust are 
treated separately. 

4.4 The evolution of an isolated turbulent molecular cloud 

Glover & Clark (2012b,c) studied the thermodynamics of turbulent 
molecular clouds with different metallicities. In particular, they per¬ 
formed calculations of 10^ M© clouds, typically modelled using 
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Figure 14. The evolution of the gas (solid), dust (dotted), and radiation (dashed) temperatures as functions of radius (left panel) and maximum density (right 
panel) during a radiation hydrodynamical calculation of the collapse of an unstable spherical 5 Mq Bonner-Ebert sphere. Each colour gives the state of the 
calculation when the maximum density reaches a different value, with the initial conditions given in black. The radiation temperature is always less than or 
equal to the gas temperature in this graph. The thick black long-dashed lines in the right panel give the evolution of the central gas and radiation temperatures 
with central density from Eig. 13 for reference (only for p > 10“^^ g cm“^). It can be seen that for a given density the gas is usually hotter latter in the 
collapse than earlier due to radiation from the protostar heating the outer parts of the core. The formation of the stellar core is apparent in the upper-left of the 
left panel (initial radius ^0.01 AU). 

SPH with 2 million particles. Their initial conditions consisted of a 
uniform-density sphere of number density nu = 300 cm“^, giving 
a radius of approximately 6 pc. The cloud was given initial ‘turbu¬ 
lent’ motions with a power spectrum P{k) oc k~^ where k is the 
wavenumber. The energy in the turbulence was initially set equal 
to the magnitude of the gravitational potential energy of the cloud, 
giving an initial root-mean-square velocity of around 3 km s“^. 

The turbulence was allowed to decay freely. The calculations were 
performed using a confining pressure of pext = 1.2 x 10^ K cm“^ 
which was implemented in the SPH equations by subtracting the 
external pressure from the pressure used in the momentum equa¬ 
tion. We use the same set up as Glover & Clark, although of course 
we are unable to use exactly the same initial velocity field. We use 
our full diffuse ISM model to obtain the results below, but we be¬ 
gin by excluding the hydrogen chemistry and heating due to the 
formation of molecular hydrogen (Section 2.5.2). In the first set of 
results, we have also excluded molecular depletion (because Glover 
& Clark did not include depletion), but since we find that switching 
depletion on or off has little effect on the results it is included in 
most of the subsequent calculations. For the dust-gas thermal cou¬ 
pling we switch to equation 37 because this is the equation Glover 
& Clark used. 

4.4.1 Results excluding hydrogen chemistry 

In Fig. 16 we plot the gas (upper panel) and dust (lower panel) tem¬ 
peratures as functions of the number density of hydrogen nuclei, 
nn, for a calculation performed at solar metallicity. The snapshot 
is taken just before the first sink particle is formed. The colours 
in the gas temperature plot indicate the average CO abundance of 
each point in temperature-density space. This is the equivalent of 
Fig. 2 in Glover & Clark (2012b) or the upper-right panel of Fig. 4 
in Glover & Clark (2012c). Since Glover & Clark did not include 
molecular depletion on to grains, we have turned this off in our cal¬ 
culation. The distribution of gas temperatures with density shows 
the same general features seen in the papers of Glover & Clark. At 
low densities (nn ^10^ cm“^) the gas temperature rises towards 



Figure 15. The gas (solid), dust (dotted), and radiation (dashed) tempera¬ 
tures as functions of radius when the collapsing spherical 5 M© Bonner- 
Ebert sphere has reached a maximum density of 0.01 g cm(lower, red 
lines), and 1.76 years later at the end of the calculation when the maximum 
density is 0.06 g cm“^ (upper, blue lines). Near the end of the calculation, 
after the stellar core forms, the radiation and dust temperatures become hot¬ 
ter than the gas temperature at radii of 400 — 3000 AU because of the 
radiation emitted by the protostar. 


~ 100 K as the number density decreases towards nn ~ 10 cm“^. 
The temperatures at densities nn ^ 100 cm“^ are primarily de¬ 
termined by the models of photoelectric heating, and cooling due 
to electron recombination and emission from ionised oxygen and 
carbon (equations 26 to 33). The gas temperatures lie beneath the 
equilibrium curve (solid black line from Fig. 1) because the equi¬ 
librium curve assumes no extinction whereas in fact the photoelec¬ 
tric heating is generally attenuated by the dust extinction (depend¬ 
ing on the location of the gas). Our maximum temperatures in this 
density regime are somewhat cooler than those obtained by Glover 
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Figure 16. Phase diagrams showing the gas (top panel) and dust (lower 
panel) temperatures as functions of the gas density. The colours in the 
gas temperature plot indicate the average CO abundance of each point in 
temperature-density space. The black solid line in the gas temperature plot 
is the equilibrium curve from Fig. 1. 


Figure 17. Phase diagram showing the gas temperature as a function of the 
gas density for a calculation that is identical to that portrayed in Fig. 16, 
except that molecular depletion is turned on. The black solid line is the 
equilibrium curve from Fig. 1. 



& Clark due to our slightly different models in this regime as al¬ 
ready demonstrated in Fig. 1. At densities nn ~ 100 — 10^ cm“^ 
there is a large dispersion in the gas temperatures with a range of 
Tg 8 — 50 K obtained at densities nu ^ 10^ — 10^ cm“^. 
In this region, although there is still signihcant photoelectric heat¬ 
ing, the low-density coolants just mentioned become ineffective. 
Instead, the main cooling is from molecular lines, in particular 
CO (e.g. Goldsmith 2001; Glover & Clark 2012b), which depends 
strongly both on density and temperature. At the same time, work 
from the gas motions becomes important (either heating or cool¬ 
ing the gas, depending on whether the gas is being compressed or 
expanding, respectively), as does heating in shocks. Together these 
effects lead to the wide range in temperatures. Finally, at densities 
riH ^ 10^ cm“^ thermal coupling of the gas to the dust becomes 
important and the gas temperature converges towards that of the 
dust as the density increases. 

As mentioned above, Glover & Clark (2012b) do not treat de¬ 
pletion of CO onto dust grains. Turning this on has little effect on 
the phase diagram (Fig. 17). As may be expected, the temperatures 
at densities of nn ~ 10^ — 10^ cm“^ are slightly higher, but the 
effect is hardly noticeable. At lower densities the CO is not de¬ 
pleted and the dominant coolant is C^ rather than CO anyway. At 


Figure 18. Phase diagram showing the gas temperature as a function of the 
gas density for a calculation that is identical to that portrayed in Fig. 17, 
except that all the carbon is assumed to be in the form of C+. The black 
solid line is the equilibrium curve from Fig. 1. 

higher densities where the CO does become depleted, the dominant 
coolant is dust 

We now explore the sensitivity of the results to several of the 
assumptions in our model. We examine the dependence on changes 
of the chemistry, in particular the abundance of C^ and the de¬ 
pletion of molecules like CO onto grains, as well as the effects of 
variations in the molecular line cooling rates and the strength of the 
coupling between the gas and the dust. 

The main effect of the carbon chemistry model described in 
Section 2.5.1 is to provide the abundances of the coolants C^ at 
low-densities and CO at high densities. Typically the transition be¬ 
tween the dominance of the two coolants occurs at nn 10“^ cm“® 
(e.g. Glover & Clark 2012b). The temperature of the gas is quite 
sensitive to this transition from C^ to CO. For example, in Fig. 18 
we give the results obtained by assuming that the carbon remains 
as C^ at all densities and we turn off all molecular cooling. In this 
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Figure 19. Phase diagram showing the gas temperature as a function of the 
gas density for a calculation that is identical to that portrayed in Fig. 17, ex¬ 
cept that the molecular line cooling rate has been increased by a factor of 3. 
The colours in the gas temperature plot indicate the CO abundance of each 
point in temperature-density space. The black solid line is the equilibrium 
curve from Fig. 1. 
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Figure 20. Phase diagram showing the gas temperature as a function of the 
gas density for a calculation that is identical to that portrayed in Fig. 17, 
except that the collisional thermal coupling between the dust and the gas is 
15 times smaller (equation 36 has been used rather than 37). Note that there 
is a larger dispersion in the gas temperatures at densities nn ^ 10^ cm“^. 
The black solid line is the equilibrium curve from Fig. 1. 

case, essentially all of the gas lies at temperatures below the equi¬ 
librium curve of Fig. 1 and there is little gas with temperatures 
greater than 20 K with densities nn ~ 10^ — 10^ cm“^. This ex¬ 
treme case shows that it is important to have at least a simple model 
that switches the carbon from to CO. 

In Fig. 19 we return to our standard model including molecu¬ 
lar depletion, but we increase the molecular line cooling rates by an 
arbitrary factor of 3. As expected, this decreases the temperatures 
of the gas with densities nu ^ 10^ — 10^ cm“^. The lowest tem¬ 
perature (at densities nu ^ 10^ cm“^) decrease from Tg 8 to 
5 K, but the highest temperatures (at densities nn ~ 10^ cm“^) 
are not significantly affected. 

As mentioned in Section 2.4.3, equations 36 and 37 for the 


gas-dust thermal coupling coefficient differ by a factor of 15 (the 
latter, used by Glover & Clark, gives stronger coupling). There¬ 
fore, in Fig. 20 we examine the effect of using the weaker coupling 
provided by equation 36 instead. As expected, this only has an ef¬ 
fect on the temperatures at high densities, nn ^ cm-3. The 

gas shows a somewhat increased spread of temperatures because 
the coupling with the dust is not strong enough to force the gas to 
adopt the dust temperature. 

In summary, the gas temperatures in solar-metallicity turbu¬ 
lent clouds are relatively independent of many of the parameters 
of our diffuse ISM model. The most important factors are the low- 
density equilibrium model (portrayed in Fig. 1) which essentially 
sets an upper limit to the temperatures at densities nu ^ 10^ cm“^, 
and the way in which the carbon transitions from C^ to CO. The 
exact molecular cooling rates, molecular depletion, and dust-gas 
thermal coupling parameters are less important. 

4.4.2 The effects of metallicity and hydrogen chemistry 

Following Glover & Clark (2012c), we examine the dependence on 
metallicity and on hydrogen chemistry. We perform full calcula¬ 
tions (i.e. including the simple carbon chemical model and molec¬ 
ular depletion) for four different metallicities: 3, 1, 1/10, and 1/100 
times the solar value. For each metallicity, we perform calculations 
that exclude hydrogen chemistry and its associated heating terms 
(results from the solar-metallicity case have already been displayed 
in Fig. 17). We also perform calculations that include the evolution 
of hydrogen (i.e. its molecular fraction xh 2 ) and heating due to H 2 
formation on dust grains. Two calculations are performed for each 
metallicity: one beginning with fully atomic hydrogen, and one be¬ 
ginning with fully molecular hydrogen. 

In Fig. 21 we plot the evolution of the fractional abundance 
of molecular hydrogen in the eight calculations that include hy¬ 
drogen chemistry up until the formation of the first sink particle 
in each calculation. This is the equivalent of Fig. 3 in Glover & 
Clark (2012c). Beginning with fully molecular hydrogen, the de¬ 
crease in H 2 is primarily due to photodissociation in the outer parts 
of the clouds. The destruction is greater at lower metallicities be¬ 
cause the UV radiation is less attenuated by dust. Beginning with 
fully atomic hydrogen, H 2 is formed on dust grains and the total 
abundance monotonically increases in each calculation. However, 
the rate of increase strongly depends on the metallicity. The main 
reason for this is that there are more dust grains on which to pro¬ 
duce H 2 at higher metallicity (equation 40) so the rate of formation 
is higher. In addition, the extra dust attenuates the UV radiation that 
destroys H 2 so the destruction rate is also lower. By the time stars 
begin to form, clouds with solar metallicity or higher are mostly 
molecular regardless of whether atomic or molecular initial con¬ 
ditions were adopted (particularly in the central regions where the 
stars actually form). However, at lower metallicities the initial con¬ 
ditions matter much more. When stars begin to form in the 1/10 Z© 
calculations, the H 2 abundance is 65% when beginning with molec¬ 
ular initial conditions, but only 15% beginning with atomic initial 
conditions. For 1/100 Z© the disparity is even stronger, with the 
two abundances being 60% and 1.7%, respectively. 

The behaviour of the molecular hydrogen evolution is quali¬ 
tatively similar to that reported by Glover & Clark (2012c). In par¬ 
ticular, the growth rate of molecular hydrogen in the atomic solar- 
metallicity calculation matches theirs to within 10%, and the de¬ 
cay rates of the abundances in all of the molecular calculations are 
very similar. The clouds in Glover & Clark (2012c) take longer to 
begin forming stars (2.7-4.0 Myr rather than 1.7-2.6 Myr). This 
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difference is most likely due to the different initial turbulent veloc¬ 
ity field. With fully molecular initial conditions, because the stars 
take longer to form in the calculations of Glover & Clark (2012c), 
the H 2 abundances when the star formation takes place are lower. 
This is most important in the lowest metallicity case, where the to¬ 
tal abundance is only 20% when the star formation begins in 
the Glover & Clark (2012c) calculation, whereas it is 60% in our 
calculation. There are also differences in the abundances of molec¬ 
ular hydrogen in the lower metallicity calculations with atomic 
initial conditions. The total H 2 abundances when beginning with 
fully atomic hydrogen remain low throughout both our calcula¬ 
tions and those of Glover & Clark. However, at any particular time 
our abundance is twice that obtained by Glover & Clark (2012c) 
with 1/10 Z© and it is about an order of magnitude higher for the 
1/100 Z© case until shortly before stars begin to form (when the 
abundances in both calculations are 1%). No doubt some of this 
difference is due to the different velocity field and dust opacities 
used in the calculations, though it is not clear whether or not this 
explains all of the discrepancy. 

Turning to the phase diagrams of gas temperature in Fig. 22, 
the general trends are again in agreement with those found by 
Glover & Clark (2012c). First we consider the calculations that ex¬ 
clude hydrogen chemistry (they assume all hydrogen is in molecu¬ 
lar form) and its associated heating (left column of Fig 22). Increas¬ 
ing the metallicity above solar results in a reduction of the lower en¬ 
velope of the gas temperatures at densities tt-h = 10^ — 10^ cm“^ 
because the additional molecular cooling and increased extinction 
produces lower gas temperatures. The upper temperature envelope 
is less affected. Lowering the metallicity to 1/10 of solar results in 
a smaller range of temperatures because extinction has less of an 
effect at low and intermediate densities so the temperatures tend 
to lie close to the equilibrium curve of Fig. 1 (which does not de¬ 
pend significantly on metallicity because the photoelectric heating 
and the cooling due to electron recombination and fine structure 
emission are all assumed to be proportional to the metallicity). At 
higher densities, the main heating sources (compressional heating 
and cosmic ray heating) are independent of the metallicity and al¬ 
though the gas is less well coupled to the dust it is still coupled 
well enough that the gas temperatures are kept around 10 K. As 
the metallicity is reduced still further to 1/100 solar, the tempera¬ 
tures move above the solar-metallicity equilibrium curve of Fig. 1. 
This is because the cosmic ray heating (which is taken to be inde¬ 
pendent of metallicity) becomes as important as the photoelectric 
heating at low metallicity. Furthermore, the magnitudes of the elec¬ 
tron recombination and fine structure cooling rates are proportional 
to metallicity, so regions that are heated by adiabatic and shock 
heating take longer to cool. 

At high densities nn ^ 10^ cm“^ where the gas is thermally 
coupled to the dust in the solar metallicity case there is not much 
change in the super-solar metallicity case. The temperatures are 
slightly cooler (about 1-2 K) because the coupling is even stronger 
and the dust has slightly cooler temperatures due to the increased 
extinction. A larger effect is seen with reduced metallicities be¬ 
cause the thermal coupling to the dust is substantially weaker. At 
1/10 solar metallicity, the gas has a maximum in the temperature 
distribution of about 15 K at nn ~ 3 x 10^ cm“^ before the dust 
cooling takes over and reduces the temperatures at higher densities. 
At 1/100 solar metallicity, the gas temperature rises with increasing 
density until nn ~ 2 x 10^ cm“^ when it reaches a maximum of 
20 K before it begins to drop again at higher densities. 

We now consider the effects of hydrogen chemistry on the gas 
temperature due to heating of the gas when molecular hydrogen is 



Figure 21. We plot the evolution of the total fractional abundance of molec¬ 
ular hydrogen in calculations of turbulent clouds up until the time that stars 
begin to form. For each metallicity, we perform two calculations: one with 
fully molecular initial conditions (upper lines), and one with fully atomic 
initial conditions (lower lines). The metallicities are 3 times solar (black 
dotted lines), solar (red solid lines), 1/10 solar (green dashed lines), and 
1/100 solar (blue dot-dashed lines). By the time stars begin to form, clouds 
with solar or super-solar metallicity are mostly molecular, regardless of 
their initial H 2 abundance, but little molecular gas has formed in the low- 
metallicity clouds that consist of fully atomic hydrogen initially. 

formed on dust grains (equation 50). This has little effect in the cal¬ 
culations that began with fully molecular hydrogen (right column 
of Fig. 22). There are two small differences. First, in the solar and 
super-solar calculations the upper envelope of gas temperatures is 
higher at densities nn ~ 10^ — 10^ cm“^. This is because some 
molecular hydrogen is destroyed by photodissociation and when it 
reforms on dust grains this adds a source of heating. Second, in the 
low-metallicity calculations some regions of mostly atomic gas be¬ 
come molecular at densities tt-h ~ 10^ — 10® cm“^ and the extra 
heating makes this gas hotter than the molecular gas at the same 
densities. Thus, a bifurcation of the gas temperature is apparent at 
these densities. 

Beginning with fully atomic gas has a much greater effect on 
the gas temperature (centre column of Fig. 22). In these cases, heat¬ 
ing due to H 2 formation on dust becomes significant at densities 
riH ^ 10® cm“®, and the gas temperature is generally much greater 
than in the calculations without hydrogen chemistry or with fully 
molecular hydrogen initially. At solar and super-solar metallicities, 
the main effect is to raise the gas temperatures in the density range 
riH ~ 10® — 10® cm“®. Above these densities the gas is mostly 
molecular and the gas is thermally well-coupled to the dust. At 
low metallicities, the gas is significantly hotter from densities of 
nn ~ 10® cm“® even up to nn ~ 10® cm“®. This is primarily be¬ 
cause of the reduced dust abundance which has two main effects: 
the atomic hydrogen persists to higher densities, and the gas does 
not become thermally coupled to the gas until much higher densi¬ 
ties. 

As mentioned above, the general temperature trends seen here 
are in agreement with those found by Glover & Clark (2012c). 
However, there are also some differences. At low-densites our tem¬ 
peratures tend to be lower because of the different equilibrium 
models that we have already mentioned (Fig. 1). But there are also 
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Figure 22. Phase diagrams of the gas temperature versus the gas density in turbulent clouds with metallicities of 3, 1, 1/10, and 1/100 times solar (top to 
bottom, respectively). The calculations in the left column were performed without any H 2 chemistry (assuming purely molecular hydrogen) and without 
heating from H 2 formation. The other two columns included H 2 chemistry and heating from H 2 formation, but the calculations in the centre column began 
with fully atomic hydrogen, and the calculations in the right column began with fully molecular hydrogen. The colour bars indicate the average CO abundance 
of each point in temperature-density space (note that the scales differ from those in Figs. 16, 17, 19, and 20). The black solid lines give the equilibrium curve 
from Fig. 1 for solar metallicity. Beginning with atomic gas results in significantly higher gas temperatures for densities nn ~ 10"^ — 10® cm“® (or even 
higher densities in the low-metallicity cases) due to heating from H 2 formation on dust grains. 
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differences at higher densities. In the solar metallicity cases, Glover 
& Clark (2012c) obtain slightly lower minimum temperatures at 
nn ~ 10^ cm“^. This may be because their molecular cooling is 
somewhat stronger (c.f. the difference that increasing the molecular 
cooling makes in Figs. 17 and 19). In the calculations at 1/10 Z©, 
the mean temperatures at each density are similar to those obtained 
by Glover & Clark (2012c) for both atomic and molecular initial 
conditions, but the scatter in the temperatures obtained by Glover 
& Clark (2012c) is greater at nn ^ 10^ cm“^. This may be largely 
due to the different density structure of the clouds (caused by the 
different initial velocity fields). In the lowest metallicity calcula¬ 
tions (1/100 Z©), the temperatures obtained by Glover & Clark 
(2012c) are generally higher than ours at nn ^ 10^ cm ^ for the 
molecular initial conditions, and above nn ^ 10® cm“® for the 
atomic initial conditions. This may be related to the higher molec¬ 
ular hydrogen abundances that were commented on when we dis¬ 
cussed Fig. 21. In our initially atomic calculation, the gas is mostly 
molecular for densities nn ^ 10® cm“® meaning that there will 
be little heating from H 2 formation above these densities. If Glover 
& Clark (2012c) have more atomic gas remaining at high densi¬ 
ties this will provide additional heating. In our initially molecular 
calculation, because the stars begin to form at 2.0 Myr instead of 
3.6 Myr, more of the initial molecular hydrogen remains. The to¬ 
tal abundance is only 20% in the calculation of Glover & Clark 
(2012c) when the temperature diagrams were plotted, whereas the 
total abundance in our calculation is 60% and all of the high-density 
gas (nn ^10® cm“®) is essentially fully molecular (i.e. there is lit¬ 
tle heating from H 2 formation). On the other hand, the temperature- 
density behaviour we obtain for the low-metallicity fully atomic 
initial conditions is in very good agreement with the results given 
in Fig. 1 of Omukai et al. (2005). For 1/10 Z© calculations, the tem¬ 
perature at high densities has a peak at nn ~ 10® cm“® at similar 
temperatures in both our calculations and those of Omukai et al. At 
1/100 Z©, both our results and those of Omukai et al. display high- 
density temperature maxima of Tg 60 K at nn ~ 10® cm“®. 


5 CONCLUSIONS 

We have presented a new method for modelling the thermal evo¬ 
lution of star-forming molecular clouds. The method combines a 
model for the thermodynamics of the diffuse interstellar medium 
with radiative transfer in the flux-limited diffusion approximation. 
The former is required to correctly model the thermal behaviour of 
molecular clouds at low densities and metallicities, while the lat¬ 
ter allows us to model protostar formation. Unlike most previously 
published star formation calculations, our new model evolves the 
temperatures of the gas, dust, and radiation separately. The code 
can also follow the evolution of atomic and molecular hydrogen. 
We have compared our method with existing literature on the ther¬ 
mal behaviour of the ISM and molecular cloud cores and generally 
obtain good agreement. 

We have also explored the sensitivity of the thermal structure 
of molecular cloud cores and turbulent clouds to many of the pa¬ 
rameters that enter our diffuse ISM model. For the gas at solar 
metallicities, the most important thermal processes at low densities 
(nH ^ 1 — 1000 cm“®) are photoelectric heating of electrons from 
dust grains and cooling due to electron recombination with small 
grains and PAHs and fine structure emission from and atomic 
oxygen. At high densities (nn ^ 10® cm“®), cosmic ray heating 
and the work done by hydrodynamical flows dominate the heating, 
while the cooling is dominated by continuum dust emission so that 


the gas and dust adopt similar temperatures. At intermediate den¬ 
sities (nn ^ 10® — 10® cm“®), most processes have a significant 
effect and there tends to be a large dispersion of temperatures in 
turbulent clouds. The abundance of changes rapidly at these 
densities and because is such an effective coolant it is impor¬ 
tant to have a model for the abundance. However, the exact 
values of the molecular cooling rates, molecular depletion, and the 
strength of the thermal coupling between the dust and the gas are 
less important. When beginning with low-density molecular clouds 
(which may not be purely molecular) the thermal evolution also de¬ 
pends significantly on relative abundances of atomic and molecular 
hydrogen due to heating from the formation of molecular hydrogen 
on dust grains. This becomes more important at lower metallicities. 

Our new method should allow more realistic radiation hydro- 
dynamical calculations of star formation to be performed, partic¬ 
ularly in molecular clouds that have low mean densities (nn ^ 
10^ cm“®) and/or sub-solar metallicities. However, we also em¬ 
phasise that this model is far from complete. In particular, it relies 
on various parameterisations of the thermal effects of many physi¬ 
cal processes, it does not calculate these processes explicitly. Fur¬ 
thermore, it contains only extremely simple models of hydrogen 
and carbon chemistries and does not explicitly model the chem¬ 
istry of other elements at all. However, the basic model developed 
in this paper could easily be extended to increase the complexity 
of the chemical modelling (at the cost of increased computational 
expense). 
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